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ABSTRACT 

A control model describing a network flow hind- 
ered by multiple pure time (or transport) delays is 
formulated. Feedbacks connect each desired output 
with a single control sector situated at the origin. 
The dynamic formulation invokes the use of differen- 
tial-difference equations. This causes the character- 
istic equation of the model to consist of transcendental 
functions instead of a common algebraic polynomial. 
A general graphical criterion is developed to evaluate 
the stability of such a problem since the literature is 
evasive in the case of multiple coupled delays. A 
digital computer simulation later confirms the valid- 
ity of such criterion. An optimal decision-making 
process with multiple delays serves as an application 
of the analytical effort. 
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SECTION I 


INTRODUCTION 

A. PROCESSES WITH MEMORY AND ANTICIPATION; THE 

DIFFERENTIAL-DIFFERENCE EQUATION 

In concept, there exists a wide class of mathematical processes which op- 
erate in the present using knowledge acquired from the past as well as informa- 
tion obtained from the future. This reveals the presence of memory, the ability 
to remember the past, as well as anticipation, the ability to predict the future. 

If these processes are to be realizable in our world then they must obey the law 
of causation: the cause of things must come before their effect. A chronological 
order is imposed and the processes, in this case, "follow their own destiny." 

The capacity to anticipate is suppressed and only present and past states are 
now considered. Despite this restriction, the realization of a man-made process 
or system* operating in the present using knowledge from the past is an elegant 
concept since it simulates the way humans learn and think. 

Processes with memory can be expressed mathematically by an interesting 
class of equations: the differential-difference equations, involving derivatives 
along with differences. The literature on such special mathematics embraces 
many disciplines and spreads over three centuries in research. The earliest 
work available is that of Bernoulli, 1 * who in 1728 was studying oscillatory mo- 
tion of strings. In applied mathematics, famous names like Lagrange 2 , Poisson 3 
and Cauchy are often associated with equations of differential- difference form. 

They are reflected by today’s contemporaries in this field such as Bellman and 
Cooke 5 , Pinney 6 , Ogutzoreli , and Wright . Phenomena possessing memory have 
also been studied in the physical sciences by Bateman 9 in radioactivity, Gerasimov 10 
in heat conduction, Lagrange 11 in sound propagation, Arley 12 in cosmic radiations, 
Minorsky 13 , Nisolle 14 , Volterra 15 , Picard 16 , and Satche 17 all in mechanics and 
elasticity. In biology, one finds differential-difference equations applied to the 
renewal process 18 ’ , natural and artificial selection , and the fight for 
survival 2 . 


*For convention purposes, appendices A and B provide basic definitions for the word system and 
its associated prefixes, such as feedback and control, in context with this work. Figure 1 illus- 
trates these definitions. 

|A superscript indicates a reference number appearing in the bibliography. 
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Economists have also made contributions with these special equations. To 
mention a few, Frish 22 in economic dynamics, Kalecki 23 in business cycles, 
James and Belz , 24 Samuelson 25 , Bateman 26 in retail trade theory, Theiss 27 
in the study of savings and investments, and Tinbergen 28 in the modeling of 
ship-building cycles are all notable names. 


B. THE DELAY PROBLEM 

In modern times, operational models with transport or time delays form 
the new vocabulary for procedures governed by differential-difference equa- 
tions. A delay means that one or more operations are temporarily suspended 
in time due to some constraints. 

In engineering, one can cite numerous references on time-delay problems: 
Tsien 29 in rocketry, also Smith 30 in electrical systems, Truxal 31 in electrical 
signals, Rogers and Connolly 32 in analog computers, Huggins 33 in systems dy- 
namics, Oetker 34 , Paul 35 , Caughanowe and Koppel 36 , and Tyner 37 all in 
process control for the chemical industry. Most recently, the time-delay prob- 
lem has found application in the theory of artificial intelligence . 38 


C. THE DELAY PROBLEM AND THE MANAGEMENT FUNCTION 

As far as management analysis is concerned, it is well-known that prac- 
tical trends in the U.S. are brought about by two schools of thought: those who 
practice the case method, relying on past cases, and those who favor the quan- 
titative approach, relying on analytical knowledge gained by studying physical 
systems. The first method is reflected by the Harvard circles while the second 
one is practiced by the Stanford and MIT groups. Forrester 39 in applying the 
latter method at MIT has stressed the importance of time delays in industrial 
dynamics such as production, inventory and sales. 

The new field of Operations Research, still part of the quantitative approach, 
renders time-delay problems imminent firstly in the network flow analysis where 
minimum transition time is usually the objective and secondly in Queueing and In- 
ventory Theory 40 . Prabhu 41 has considered similar mathematics in his work on 
models for dams and grain storage. 


D. THE GENERAL OBJECTIVE OF THIS RESEARCH 

In the real world, autonomous organizations are confronted with the daily 
task of applying continuous improvement mechanisms in order to enhance the 
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viability of their working processes. Such strategy involves tactical efforts on 
the part of decision-makers. Faced with commitments and schedules, the 
management firstly evaluates the input resources at hand, secondly applies a 
planned intervention, and finally seeks feed-back, which is often delayed in time, 
in order to determine the effectiveness of the motion taken. 

Such a philosophical management model has possible application in the 
following: 

1. Transportation networks suffering delays in schedule. 

2. The administration of federal, state or local programs where a compli- 
cated bureaucracy often delays needed resources such as monetary 
funds. 

3. Administrative control of the economy such as the imposition of ceilings 
on prices and wages. In this case, feedback of the results created is 
heavily delayed due to the complexity of the economic sector. 

This paper intends to highlight the above applications through mathematical 
modeling. The effort illustrates the decision-making process involved in a 
network flow containing information and resource channels hindered by coupled 
pure time delays. 
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SECTION II x 

MATHEMATICAL BACKGROUND 

A. ASSUMPTIONS LEADING TO THE DIFFERENTIAL-DIFFERENCE 
EQUATION 

One can interpret processes with memory and anticipation using the follow- 
ing general functional equation: 


® [y (x, t)] = f {y (t), y (t ± tj (t)), 


y (t ± t 2 (t)), - • y (t ± t n (t)), 
x (t), x (t ± tj (t)), x (t + t 2 (t)), ... 


( 1 ) 


* (t ± t n (t))>, 

where $ is a generalized n th order (n > 0) Jlnear or non-linear, differential or 
difference, partial or ordinary operator, y and x are vectors of the various 
dependent and independent variables (functions of time) respectively, tj, t 2 , 

. . . t represent backward or forward delays in time, themselves possible 
functions of time. Terms like y (t ± t n (t)) signify that y is evaluated at a time 
t shifted in the future by +t n (t), or in the past by -t n (t). Such a formulation 
must also be accompanied with proper initial, boundary and constraint conditions. 


Unfortunately, no tools are available to fully treat eq. (1). Simplifying as- 
sumptions must be made. If one is restricted to a realizable, i.e. non-anticipatory, 
linear, i.e. the principle of superposition holds, deterministic (non-probabilistic) 
system and furthermore, if the dependent variables are functions of time only, 
and having constant coefficients along with constant delays then eq. (1) reduces 
to the following differential-difference equation: 


D [x (t)] = f {x (t - tj),x (t - t 2 ), • • ■ x (t - t n ),x (t), t> (2) 

Preceding page blank 
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where 


D [ ] 



d n d n ” 1 

+ a i — + 

n-1 1 x 

dt n dt"- 1 


+ a i 


_d_ 

dt 


is the linear ordinary differential operator, with the a’s being constant coefficients. 

If, in addition, the process is instantaneously dependent on just the present 
state, then eq. (2) further reduces to the ordinary differential equation which 
traditionally governs countless mathematical models: 


D [x (t)] = f (x, t). 


(3) 


Eq. (3) is indeed restricted when compared to eq. (1). 

Equation (2), which assumes present and past states, is expressed explicitly 
for a single input-output system as: 


d" 

a n x (t) + a 

dt n 


in-l 


n-1 




x(t) + ••• + a 


l ^ x (0 = b o x ( t )+b 1 x(t-t 1 ) + b 2 x(t-t 2 ) 

(4) 


+ * • * + b n X ( t ~ t n'> + U (t). 


where u(t) is a prescribed input activating the process. The a's and b's are 
constant coefficients. 


B. THE LAPLACE TRANSFORM OF A DIFFERENTIAL-DIFFERENCE 
EQUATION 

Since Eq. (4) is linear it can be transformed to a simple polynomial form 
using the well-known Laplace Transform, £ . Remembering that 


£ 


df(t) 

dt 


s and £ [f (t - tj)] = e 


6 


C 



Eq. (4) becomes: 


+ a 


L~1 


+ a . 


s = b n 


-sti + b e“ st 2 


b e 


-st r 


+ £{u(t)} ( 5 ) 

(e = 2.7- • •) 


One obtains the characteristic equation of the system that eq. (4) governs by 
rearranging eq. (5) and suppressing the input. Therefore, 


F(s, e s ) = a n s n + a n _j s 


,n _ 1 


s - b 0 - bj e -st 1 - b 2 e -st2 


b e _st n = 0(6) 


is called characteristic to the system because it describes the free behavior, 
the outside input being suppressed. This equation helps determine stability, a 
salient feature in any system study, as later shown. 


C. THE MATHEMATICAL FORM OF THE PURE TIME DELAY 


Terms containing e -s in eq. (6) correspond to the various transport lags or 
delay terms such as y (t - t £ ) present in eq. (4). The syntax of the word "time- 
delay” is often misleading due to the existence of numerous synonyms. Depend- 
ing on the field of study, the following definitions occur in the literature: 

(1) distance-velocity lag (Mechanics) 

(2) retarded action (Mechanics) 

(3) dead time (in Electrical Engineering) 

(4) transport lag or delay (in Transportation Science) 

(5) pipe-line delay (in Chemical Engineering) 

(6) iddle period (Operations Research) 
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(7) pure time delay (Applied Mathematics) 

(8) retarded argument (Applied Mathematics) 

(9) translation operator (Applied Mathematics) 

They all refer to the same transformation: 

pure-time delay 

y (t)— ► -*-y (t - t t ) 

, of period tj 

The confusion arises when another function is mentioned: the phase or time lag. 
The latter is called, at a loss, an exponential delay. 39 Various such transforma- 
tions are illustrated in Figure 2. The first-order time lag l/(sr+ 1), for in- 
stance, produces an output at the instant the input activates it. The pure time 
delay e _stl , on the other hand, does not respond for a time period tj after the 
input activates it, and only then, the response comes out as the exact image of 
the input. 

This work will adopt the following notation: 


— S t 

£[f(t-t i )]=e 1 = pure time delay (which concerns 

most of this research) 

— ^ / fj- 1 

£ [e *] = = simple time lag (1 st order) 

s r 1 + 1 
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SECTION III 


THE NETWORK FLOW GRAPH METHOD FOR 
PROCESSES WITH PURE TIME DELAYS 

A. FROM THE DIFFERENTIAL-DIFFERENCE EQUATION 
TO A FLOW GRAPH 


Instead of working directly with a differential-difference equation it is pos- 
sible to manipulate its Laplace equivalent, in a graphical form. As an illus- 
tration, consider the following coupled system of differential-difference equations 
having a pure time delay and a simple time lag: 


a prescribed function, f (t), 
serving as an i npu t . 

(7) 

*2 (O = *i (0 - 3 x 6 (t) 

(8) 

*3 (t) = K x 2 (t) 

(9) 

*4 (t) = e" t/T x 3 (t) 
(time lag) 

(10) 

(t) = x 4 (t) = j* x 4 (t) d t 

(11) 

: 6 (t) = X 5 ^ * t l> 

(pure time delay) 

(12) 


(°) si gnifi es d/dt. K, A, r and t x are constants. The Laplace transforms to 
equations (7-12) are: 


X x (s) = F (s) 


(13) 
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s x 2 (s) = s X, (s) - a x 6 (s) 


(14) 


S X 3 (s) = K S x 2 (S) 

(15) 

s (s) = 1 • s X, (s) 

(16) 

" r s + 1 

X 5 (S) = x 4 (s) =!• s X 4 (s) 

(17) 

“St, 

x 6 (s) = e ‘ x 5 (s) 

(18) 


Utilizing the flow graph method exposed in Appendix C, eqs. (13-18) can be 
represented by the network flow appearing in Figure 3a. Such a graph consists 
of three major terms, shown in a brief version in Figure 3b as: 


H (s) = H 


K . 1 

r s + 1 s 


(the series product of all 
<terms in cascade excluding 
|the pure time delay, 



(19) 


D (s) = D = e 


the pure time delay, 


( 20 ) 


G(s) = G=-a=the feedback loop. 


( 21 ) 


B. THE TRANSFER FUNCTION FOR A SINGLE PURE 
TIME DELAY PROBLEM 

Applying Mason's Loop Rule (see Appendix C) it is possible to obtain the 
transfer function (T.F.) defined in Appendix B for the system in Figure 3b as 
follows : 
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Pj = H D 


(22) 


L t = - H D G (23) 

A=1-L x =1+HDG (24) 

Aj = 1 (25) 


where P t is the only forward path, L x is the only feedback loop factor, A is the 
determinant, and A is the cofactor. Thus 


T. F. 



l 


Pi A t 

A 


H D 

1 + H D G 


(26) 


Setting the denominator of equation (26) equal to zero, yields the characteristic 
equation for the process: 


1 + H D G = 1 + H (s) -D (s)*G (s) 


= 1 - 


a K 


1 

• — • e 


(rs + 1) s 


1 = 0 


(27) 


(C) THE CASE OF MULTIPLE TIME DELAYS: 

Consider a more complicated case. One can allow for two pure time delays 
as shown in Figure 4a. For this case one has, again for Mason's Loop Rule: 


Pi = H, D, 

(28) 

P 2 = h 2 d 2 

(29) 

Li = - Hj Dj G x 

(30) 
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L 2 - H 2 d 2 g 2 


(31) 


A = 1 - (Lj + L 2 ) = 1 + Hj Dj Gj + H 2 D 2 G 2 


(32) 


A 


i = 1 


(33) 


A 


= 1 


(34) 


and the T.F. is: 


T. F. 


Hj H2 D! D 2 

1 + H x Dj Gj + H 2 D 2 G 2 


(35) 


This problem can be said to be a trivial extension of the single delay case. This 
is so because it can be separated into two single delay problems as shown in 
Figure 4b. The two loops are not coupled since the output of the first system is 
directly used as the input to the second one and the latter has no feedback in- 
fluence on the first one. If coupling is present however then the system is not 
trivial. Figure 4c shows the non-trivial case for two delays. For this flow 
graph, one has: 


Pj = H x H 2 D x D 2 (36) 

P 2 = 0 (37) 

L i = - H i D i G i (38) 

L, = - Hj H, Dj D 2 G 2 (39) 

A = 1 - (L j + L 2 ) = 1 + Hj D 1 (Gi + H 2 G 2 D 2 ) ( 40 ) 

A 1 = 1 (41) 
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and therefore, 


zx*. 


T. F. 


Hi H 2 Dj D 2 


1 + H x D 1 (H 2 G 2 D 2 + G x ) 


Extending the work to a non- trivial problem with N delays, the flow is 
shown in Figure 4d. Its general T.F. can be shown to take the form: 


IN 

TT D i 


T. F. = 


1 + £ hit H i 

j=l L i=l 


where 


L = summation of terms over the j index from 1 to N, 
j = i 
N 

fj = product of terms over the i index from 1 to N. 

i=l 

— s t . 

Dj = e 1 (pure time delay) 

G. = feedback terms 

H £ = product of all terms in each i th forward path between nodes, ex- 
cluding each D. . 

If one sets the denominator of equation (43) equal to zero then the result yields 
the characteristic equation for the general non-trivial case: 


F (s, e s ) = 


N i- j -i 

= 1 * J2 G| TT h ‘ d ‘ = ° 

j = 1 _ * i= 1 


(43a) 


with G . and Fi containing algebraic terms in s, and Di containing the exponential 
terms e” sti . As an illustration, if three non-trivial pure time delays are 
present then equation (43a) takes the form for N = 3: 


F (s, e s ) = 1 + Gj Hi D x + G 2 Hj H 2 Dj D 2 + G 3 U l H 2 H 3 D 2 D 3 


= 1 + G x (s) H x (s) • e 1 + G 2 (s) H x (s) H 2 (s) • e 


•Sl(t,+t 2 ) 


+ G 3 (s) Hj (s) H 2 (s) H 3 ( s) • e 


-s(t 1 + t 2 + t 3 ) 


(43b) 
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SECTION IV 

THE STABILITY OF PROCESSES 
HAVING PURE TIME DELAYS 


A. DEFINITIONS 

A linear system is defined as stable if the time function representing its 
response to a simple impulse input remains bounded in amplitude as t — 00 , and 
unstable otherwise. More important, a stable system is asymptotically stable 
if its impulse response tends to zero as t — 00 . Since the transient response 
resulting from an arbitrary input, or disturbance, can contain only time func- 
tions of the type which occur in the impulse function, it is a corollary to this 
definition that an asymptotically stable linear process will eventually return 
to equilibrium for any transient disturbance. It can also be shown that the out- 
put of an asymptotically stable process will remain bounded if the input is 
bounded. In the real world, processes are not created to simply be stable for 
a short period. They must be stable in the long run. Therefore asymptotic 
stability (t - ®) is the prime objective. To shorten notation, in the rest of this 
work the word stable will signify asymptotically stable. Finally if a system is 
stable one desires to know how close it is to being unstable. This is the addi- 
tional concept of relative stability as compared to absolute stability, previously 
discussed. 


B. SURVEY OF THE LITERATURE ON STABILITY: 

Numerous criteria based on the previous definitions are available to 
provide an answer to the question of stability for linear systems. These methods 
are either analytical or graphical. Basically, they all start with the character- 
istic equation derived by Laplace transformation. Traditionally, the analytical 
approach is reflected in the Routh and Hurwitz criteria and the Root- Locus 
method. The first two answer the question of absolute stability while the last 
one determines relative stability. On the other hand, names like Nyquist, and 
Bode often appear when the graphical approach is selected. Numerous other 
criteria have also been derived. 42 


Preceding page blank 


15 



(i) The Analytical Methods: 


The Routh and Hurwitz criteria give fast results 43 for the following 
characteristic equation: 


a Q + aj s + a 2 s 2 + a 3 s 3 + . . . a n s n = 0 


( 44 ) 


where the a.’s are constants. This equation corresponds to a process governed 
by ordinary differential equations with no differences, and where n, the number of 
roots of equation (44), is finite. 

The exponential terms in equation (43a), on the other hand, are periodic by 
definition: 


"stj “(cr + icd)t 1 t j " i CO t j 

e = e - e * e 


- i (co t J + 2 n7T ) 


(n = 0, ±1, ±2. . .) (45) 


(i = >TT"I) 


and thus equation (43a) is transcendental with an infinite number of roots. 

Although the stability requirements are the same for a purely algebraic 
characteristic equation as well as for one with exponential terms, the finite 
analytical criteria fail to give exact results for the latter. This is due to the 

presence of the infinite order polynomial which represents the exponential as 
follows: 


~ S 

e 


= 1 - tj s + 


t 2 s 2 
2 ! 


3 ! 


(46) 


Attempts have been made to approximate e stl , by truncating this Taylor series 
expansion. This approximation has been shown to be contradictory. Choksy 44 
demonstrated that, depending on where one truncates the series, the process 
can be found to be both stable and unstable when the Routh criterion is later ap- 
plied, for instance. Other approximations to e _st >, have been derived in order 
that the Routh criterion can handle it. Truxal 31 discusses various ones 
in increasing order of sophistication. Firstly, he starts with the exponential 
function expressed as a limit: 
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( 47 ) 



but mentions its weakness at n = 3. Secondly, he proposes the McLaurin series 
approximation as: 


t? s 2 


1 + t x s + 


2 ! 


t 3 s 3 
3 ! 


(48) 


which is derived from the previous Taylor series expansion, again with limi- 
tations. Thirdly, the Pade approximation table is presented. This is a ratio of 
two finite polynomials with selected coefficients (see Figure 5). One can choose 
any one fraction in this table as an approximation to e~ s . The higher order poly- 
nomials in the numerators or denominators produce a better accuracy. The 
Pade approximation and its extension, the technique of Single and Stubbs 
which produces a similar ratio, are favorite methods in industry. ’ 


Instead of customizing the troublesome exponential to the^Routh criterion, 
an effort has been directed to the root locus technique. Krall ' was initially 
interested in the roots of transcendental equations from an involved mathematical 
point of view. Finally in 1967, he produced a digital algorithm 48 used to generate 
the root locus diagram for an equation having an e' s term. Such studies were 
pioneered by Chu 49 in 1952. 

All of these works, so far, have considered only the case of the single pure 
time delay, thus providing restricted results. Contributions to the problem of 
multiple delays are very rare, probably due to the fact that most real-world 
applications assume that the multiple case is a trivial extension of the single 
case. This is not always so, as discussed in Section IE. Yuan-Yun, Quing and 
Lian 5 0 produced a paper on multiple delays from a sophisticated analytical 
view. In addition, Shaughnessy and Kashiwagi 51 have derived a stability indica- 
tive function for the multiple delay case. 


(ii) The Graphical Approaches: 

These methods, such as the Nyquist criterion, can determine exactly the 
relative stability of systems with pure time delays. They handle the latter 
without any numerical approximations. The Nyquist criterion and its application 
to the case of a single delay are outlined in Appendix D. Figures 6a and 6b refer 
to the text in this appendix. 


17 



C. THE SATCHE TECHNIQUE: 

It can be seen that the simple example used in Appendix D produces a Nyquist 
plot difficult to draw. The diagram is even more troublesome to evaluate for 
higher order polynomials, not to mention the multiple non-trivial cases. 

Satche 1 7 derived the following elegant procedure which greatly simplifies the 
Nyquist plot. One supposes that the characteristic polynomial in equation (43a) 
can be separated into pure algebraic and pure exponential parts. Therefore, 

F(s, e s ) can be written as: 


F (s, e s ) = F x (s) - F 2 (e s ) 


(49) 


The loci for F x and F 2 are separately drawn in Figure 6c according to the Ny- 
quist method. Looking at that figure, the Fj (s) curve is a closed form algebraic 
polynomial in s. It closes on itself and part of it resides at an infinite distance 
from the origin. The area it encloses is cross-hatched. The F 2 (e s ) locus, on 
the other hand, describes a finite contour around the origin. The area enclosed 
by F 2 is covered with dots. In the Nyquist Diagram, one recalls that it is 
necessary to determine the number of rotations performed by a vector, originat- 
ing at the origin, whose head follows the contour F as shown in figure 6a. This 
number of revolutions is identical to the number of roots with positive real 
parts (unstable roots). In the Satche diagram (Figure 6c), a similar vector is 
drawn for each contour in F x and F 2 . Then, according to equation (49), F(s, e s ) 
is the resultant vector for F t - F 2 . It joins the points A and B in the diagram. 
Both ends of the vector AB are now free to move. The reader can follow the 
motion of AB as OA and OB are rotated clockwise in Figure 6c. AB revolves 
on itself while its magnitude changes. At the lower intersection of Fj and F 2 
the magnitude is zero and as soon as the intersection is passed, AB emerges 
with a complete half rotation. The same phenomenon occurs at the other inter- 
section. Thus the intersection points are critical. Each point forces AB to 
rotate half a revolution. The Nyquist criterion can now be modified to impose 
on AB the condition that if the vector performs one or more complete revolution 
on itself then the system is unstable. In summary, three cases arise in this new 
criterion called the Satche criterion: 

(1) If the two curves in the Satche diagram are completely disjoint, the 
system is stable. This case is obvious since the angle of the rotation 
of AB is bounded and always less than 360°. 

(2) If one curve completely encircles the other without intersection, the 
system is unstable; this is another obvious case, since at least one 
revolution of AB is inevitable. 
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(3) If the two curves intersect, then stability can be determined as follows. 
Recalling that s = i" in the Nyquist plot, let ^ a) - oj 2 be the param- 
eter range of two successive outermost intersection points on the 
curve of Fi during which Fi lies within the area bounded by the 
F 2 curve, (i.e. the area where instability can occur corresponding 
conformally to the right-half of the complex s-plane in Figure 6b.) 

Then if that part of the curve of F 2 corresponding to the same range 
lies completely outside the area bounded by F , then the vector AB 
suffers no net rotation and the system is stable. 

This case is demonstrated by Figure 6c as follows: One notes that the 
vector A 2 *B 2 just above the first intersection proceeds to become A 3 B 3 
without any rotation, thus indicating stability. Once the upper inter- 
section is crossed, A 3 B 3 rotates 180° to become A 4 B 4 . The latter 
proceeds clockwise, and somewhere after AB, it eventually has the 
same direction as A^B 3 indicating one complete revolution. Any 
further rotation dictates instability. 

To the author's knowledge the above Satche criterion has been applied only 
to single pure time delay cases. It is possible to extend the criterion for general 
characteristic equations with numerous exponentials provided that the algebraic 
and exponential terms can be separated. The proof of this extension involves a 
simple axiom from vector algebra, namely the law of uniqueness in vector ad- 
dition. This states that to every pair of vectors, there is a unique vector, called 
their sum. As an illustration, equation (49) can be written in detail as: 

i 

F (s, e s ) = Fj (s) - F. ( e" s 1 i - 1 ) ( 50 ) 

j = 2 

where F 2 , F 3 , . . . F i form a series of curves on the Satche diagram, each one 
corresponding to a different exponential in equation (43a). One can add up the 
vector joining a point on the F. locus to a corresponding one on the F._ locus. 
The result is a single vector (the addition law) which corresponds to a new locus, 
F i + F A _ . . The new vector is then added to the appropriate vector for F i _ 2 . 

The procedure is repeated until a final curve and corresponding vector are ob- 
tained con taining all of the exponentials. This curve is then compared to F x (s) 
as done in the simplified Satche diagram, in Figure 6c. For the case of numer- 
ous exponentials, the F 2 curve has more than one intersection with F t . For 
stability, one should consider only the outermost crossovers. 
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D. COMPUTER IMPLEMENTATION 


The Generalized Satche criterion is readily adaptable to computer pro- 
gramming. Such effort is performed in order that the computer can determine 
numerically whether a system is stable or not without recourse to a diagram. 
The Satche diagram thus becomes secondary and serves only as a graphical 
summary of the numerical work. 

Proper instructions in FORTRAN are set up for an IBM 360/91 system. 

The program, accepts any number of pure time delays, within the capability of 
the machine. The graphical capabilities of CALCOMP (Computer Graphics) are 
then used and the final result is a Satche diagram for a non-trivial multiple 
pure time delay problem drawn by the machine. The application that follows 
illustrates the procedure. 
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SECTION V 


APPLICATIONS 


A. AN OPERATIONAL MODEL WITH TWO NON- TRIVIAL PURE 

TIME DELAYS: 

Consider a real network possessing both an information and a material, or 
resource, flow. Information flow is almost instantaneous since it usually in- 
volves electrical or electronic means such as telecommunications. Resource 
flow, on the other hand, suffers delays in schedule such as those in transpor- 
tation and in office bureaucracy. The process is shown in Figure 7 and con- 
sists of: 

1. A management sector made up of: 

(a) A single junction point where feedback information from two desti- 
nation points is subtracted from the forward resource flow rate 
entering, (data processing), 

(b) A control having the decision-making option to vary a factor, K, 
which, in turn , takes on fractional values between 0 and 1 and 
multiplies the resource flow rate such as to increase or diminish 
it. 

2. An undesirable forward time lag (1 /Ctj s + 1] ) which deforms the char- 
acteristics of the resource rate function. The time-lag occurs quite 

39 

often in the real world, as demonstrated by Forrester. 

3. An integrator which changes the resource flow rate into a flow of 
individual units. 

4. Two pure time (or transport) delays, hindering the flow which eventually 
reaches two destination points. 

5. An information line connecting each of the two destinations with the 
original junction, bringing back knowledge on the volume level in time- 
past due to the presence of the delays. 

The input takes the form of a rate (units/time) since the m ana gement sector 
traditionally deals with rates. The destination points on the other hand, situated 
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at a lower level in the organizational hierarchy, usually report on the volume 
level of actual units received. An integration is therefore required to transform 
the rate of flow into a level of flow. Figure 8 represents the mathematical flow 
graph for Figure 7. One notes that the model constitutes a non- trivial case since 
the two destination points report to the same centralized data processor. The 
differential-difference equations extracted from Figure 8 are: 


4 - b l X 6 " b 2 X 7 CO 

(51) 

x 3 <*) =Ki 2 (t) 

(52) 

* _t / T 1 * 
x 4 (*) = e x 3 (t) 

(53) 

x 5 (t) = j x 4 (t) d t 

(54) 

X 6 (O = X 5 ( t ~ fc l) 

(55) 

X 7 <*) = X 6 ( 4 ‘ t 2 ) 

(56) 


where 



*! ( t) = an arbitrary input rate (units/time) originating from the environ- 
ment. 

(t) - the resultant rate exiting from the junction point (units /time) 

x 3 (t) = the rate of resource flow (units/time) resulting from the manage - 

ment control on x 2 (t). 

*4 (*) = the rate of resource flow (unit/time) transformed by the undesirable 

time-lag. 
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x g (t) = the level of resource flow (units), once the above rate (units/time) 
is integrated. 

x ( t ) , x (t)= the levels (units) at the two destinations during which the flow 
suffers pure time delays of tj and t 2 , respectively. 

K = the management control (dimensionless) which is to be varied. 

b x , b 2 = the feedback constants (1/time). 

Equation (52) describes the transformation at the central junction; feedback 
terms such as b x x g (t) and b 2 x ? (t) (units/time) represent the rate at which 
the resource units entering the destination points are being rejected due to 
overflow created by constraints. A linear relationship has been assumed be- 
tween the rejection rates b l x g (t) and b 2 x ? (t) and the levels x g (t) and x ? (t), 
since a higher "crowding" of units at these receiving stations causes the latter 
to produce a higher rejection rate. The central junction processor (such as a 
computer) is programmed to take the above into account and thus pre-plan the 
schedule by subtracting the rejection rates from the entering resource rate as 
shown in equation (52). That difference enters the management control where 
an optimal value of K has to be determined. The latter should allow for maxi- 
mum entering flow rate and, at the same time, provide asymptotically stable 
level fluctuations at the destination points. 


B. APPLICATION OF THE GENERALIZED SATCHE TECHNIQUE: 

To determine such optimum value of K, the generalized Satche technique is 
used. Apply ing first the Laplace Transform to equations (51-56), six equations 
are obtained in the s domain. Combining these into one and rearranging terms 
as in equation (5), the characteristic equation becomes: 


F (s, e s ) = s (r, s t 1) + bj K e Stl + b 2 K 


’(tj+tj) 


(57) 


The same result can be obtained directly from the graph in Figure 8 along 
with the general equation (43a) for N = 2. The latter becomes: 


F (s, e s ) = 1 + Gj H x D x + G 2 H x T> 1 H 2 D 2 


(58) 


= 0 
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and since in Figure 8, 


G i = - b i 


G 2 = - b 2 


H, = K 


Tj S + 1 S 


H, = un i ty 


“ S t i 

Dj = e 1 


_ "St. 

D 2 = e 


then equation (58) eventually becomes equation (57). Since equation (57) takes 
the form of equation (50), one has: 


F (s, e s ) = F x - F 2 


(59) 


with 


and 


Fi = s (Tj s + 1) 


F 2 = - bj K e 1 - b 2 K e 


s ( t l +t 2) 


(60) 

(61) 


Subst ituti ng s = i a; (see Appendix D) into equations (60) and (61) one ob tains 
(i = /^l): 


Fj (O)) = - Tj OJ 2 + i 


CO 


(62) 
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F, (co) = 


- bj K e 1 - b 2 K e 1 2 ' 


= - bj K {cos (co tj) - i sin (co tj)} 


- b 2 K {cos [( t j + 1 2 ) - oj] — i sin [(tj + t 2 ) • co] } (63) 

Separating real and imaginary parts yields : 

X! (co) = Re Fj (co) = - Tj 0)2 ( 64 ) 

Y x (co) = Im F x (co) = co ( 65 ) 

X 2 (<^>) = He F 2 (co) = - bj K cos (co tj) - b 2 K cos [(tj + 1 2 ).-o>] (66) 

Y 2 (^-') = Im Fj'(w) = bj K s in (<d t t ) + b 2 K s in [( tj + 1 2 ") . co] 

where Re = "real part of" and Im = "imaginary part of." The constants in 
equations (64-67) are arbitrarily selected as: 


b l = b 2 = 1 

(68, 69) 

ti = 1 

(70) 

b 2 = 2 

(71) 

Tj = 1.0 

(72) 

K = 0.45 

(73) 


remembering that 0 £ K £ 1. 
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c. COMPUTER APPLICATION OF SECTION (B) 

The general program described in Appendix E is utilized to obtain the neces- 
sary numerical and graphical results for stability. The procedure consists of 
the following five major steps. 

(i) Generation of Values for Xj, Y x , X 2 , Y 2 in Equations (64-67): 

The parameter co is varied over an interval of to +rr with an incre- 
mental value arbitrarily chosen as -n / 256. Thus co takes on 513 values. Cor- 
responding numerical values for Xj, Y p Y 2 , and Y 2 are generated by the com- 
puter utilizing equations (64-67). The output is shown in Figure 9. There, the 
top headings state the numerical values of the constants as in equations (68-73). 
"B(l) TO B(M)" and "T(l) TO T(M)" represents equations (68-71). Subscripts 
appear in parantheses, with M = 2, the number of pure time delays in this case. 
Due to the presence of sine and cosine in equations (66-67), X 2 and Y 2 are 
periodic and thus redundant outside the + v to -v interval inw. The increment 
for co is selected such as to produce satisfactory smooth curves when computer 
graphics are later used. 

The nature of equations (64-65) dictates that Fj is a parabola symmetric 
about x-axis, open to the left and going through the origin on the Satche diagram. 
It closes on itself at infinity as in Figure 6c. The area of the right of this parab- 
ola is enclosed by F t . Equations (66-67) indicate that F 2 is a closed curve 
around the origin and thus Fj and F 2 intersect once or more in the Satche 
diagram. Therefore, of the three possible cases outlined in the Satche criterion 
in Section IV, only the last and most general one will be examined here. 

(ii) Locating the Regions of Intersection of F i and F 2 : 

The set X 1 and Yj along with the set X 2 and Y 2 each corresponds to an 
ordered pair (a point) on the F, and F 2 curves, respectively. It is first neces- 
sary to determine numerically the points where these two curves intersect. 

The SCAN SUBROUTINE, a numerical algorithm (see Appendix E), performs this 
required step as follows. Starting with the first line of X^ and Y in Figure 9, 
the entire ordered X 2 , Y 2 set is scanned for that pair on F 2 closest in distance 
to that point (Xj, Y } ) on Fj. The algorithm utilizes the distance definition for 
two ordered pairs: 


D = VCXj - X 2 ) 2 + (Y[ - Y 2 ) 2 (74) 


26 



The pair X 2 , Y 2 found closest is at a distance identified as DMIN. Further- 
more, if DMIN is less than or equal to the value 0.025 (an estimated very small 
term) then X 2 , Y 2 is not only the closest to X x , Yj but at these four values 
DMIN is near zero and F x and F 2 intersect. The next pair (X , Y x ) on F in 
Figure 9 is selected and the above procedure is repeated until all of the points 
on F ! have been examined. Figure 10 gives the results of the work done by 
such algorithm. The code N = 0 signifies that X 2 = ReF 2 and Y 2 = Im F 2 , on that 
line, form the closest non-intersection point resulting from the scanning of the 
entire set of X 2 , Y 2 . Furthermore, N = 1 indicates a region of intersection. One 
notes that Fj and F 2 do not have each the same co at N = 0 or 1. The examination 
of the column for N in Figure 10 reveals that a long series of 0's first appear 
followed by a few l's, and then more 0's. These l's indicate a region of inter- 
section, with probably the middle line, within these l's, as being the desired 
point. The fact that a region and not a point of intersection can be determined 
is explained by the small value 0.025 taken as criterion. The smaller it is, the 
smaller the l's region becomes and the better the resolution. However, an 
estimated value that is too small, such as less than 0.025, can cause the scanning 
algorithm to miss some intersections depending on how close to orthogonal 
(perpendicular) the meeting of F t is to F 2 . 

(iii) Locating the Outermost Intersections: 


The computer program next examines which two regions of intersections 
are the outermost ones with an increasing co for F . These two points are 
identified in Figure 11. The last column in that figure gives the bounds for co 
corresponding to F 2 . The numerical values under "OMEGA OF Fl" in that figure 
correspond to the lower and upper bound on co for F 2 and refer to co t and &> 2 , 
respectively, in the Satche criterion. 

(iv) Determining if the Values of F 2 (coj through F 2 (co^ are Enclosed by Fj : 


The set of pairs X 2 (o>j), Y 2 (o>j) through X 2 (co 2 ), Y 2 (co 2 ) represent the 
critical part of F 2 and must be examined for possible intersection with Fj . 
Subroutine SCAN is again utilized, but this time F x and F 2 are replacing each 
other in the algorithm, namely, the entire F } is scanned for each selected point 
X 2 , Y 2 within F 2 (ajj) to F 2 (a > 2 ). If an intersection occurs then a section of that 
critical part of F 2 must be to the right of F and instability is dictated. If no 
intersection occurs then two possible cases arise, either F (co t ) through F (a> 2 ) 
is completely to the right or completely to the left of F t . If it is to the right, 
the system is unstable, if to the left, stability is declared, according to the 
Satche criterion. The computer program performs the above logic and pro- 
duces output of this step in Figure 12. One should read the explanatory N and 
L codes in the figure. 
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(v) The Satche Diagram: 


The program next stores the preceding results (Figures 9 through 12) 
on a magnetic tape. Subroutine SATCHE performs this step such that the tape 
can be later mounted on a Calcomp (California Computing) System. This hard- 
ware electronically deciphers the code and mechanically draws the Satche dia- 
gram. The latter is shown in Figure 13 and summarizes the complete effort 
for the case of K = 0.45. One can observe five intersections between the para- 
bola for Fj , identified with ++++, and F 2 . F t closes on itself at infinity in the 
far right plane. That part is not drawn since it is of secondary importance. 

The critical part of F 2 (oj 1 ) through F 2 (^ 2 ), identified by **** goes to the right 
beyond the outermost intersections and therefore the case is unstable. Figure 
14 shows a highly unstable case with K = 1.0. In this instance the resource flow 
rate is at a maximum with no management restraint. Figure 15 indicates a stable 
system for K = 0.25. Figure 16 shows a case where for K = 0.35 the critical 
part of F 2 (Wj) through f 2 K) is to the left of F x and terminates at the outer 
intersection points. The value K = 0.35 is therefore an optimal case since it is 
the largest K for which the process is still stable. 

It is possible to confirm the integrity of the generalized Satche criterion » the 
computer program and the plotting routine by performing a relatively fast digital 
simulation. Since this approach falls outside of the stability theory presented 
so far it becomes a fairly rigorous test. The Continuous System Modeling Pro- 
gram (CSMP) package provided by IBM 52 , performs such simulation. Figure 17 
shows the CSMP model which corresponds to the flow graph in Figure 8 . The simulation 
produces the time response of the process when an input is applied. It solves the 
system of differential-difference equations (51-56) utilizing a digital integration 
method, the Adams technique, provides a discretized evaluation of the two pure 
time delays, and finally produces a numerical and graphical display of the out- 
put to the system. In this model simulation, a unit step function is arbitrarily 
chosen as input. One is interested in the fluctuations as time progresses of x ? , 
the final destination level, representing the output. Figure 18 shows the results 
for K = 0.45. First of all the response of x ? does not appear until time = 3.00. 

This represents the overall initial pure time delay of t x (= 1) + 1 2 (= 2) = 3. 
Furthermore the level at x ? in Figure 18 shows growing oscillations indicating 
an unbounded output thus causing the overall process to be unstable. Figure 19 
shows a stable condition when the above simulation is repeated with K = 0.25, 
this time. Here, the final destination level fluctuates significantly at the start 
but such transient behavior eventually decays to a settled level as dictated by 
the control. The stability for this process is evident at the end of 80 time units. 

The simulation is performed a final time with K = 0.35. The oscillations for 
this case are shown in Figure 20. The oscillation heights slighly decrease 
from one peak to the other revealing eventual asymptotic stability but with a 
"settling" time that is much longer than for K = 0.25. The decision-making con- 
trol is confronted with two alternatives. Depending on the trade-off objective, 
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the management can apply a tight control on the resource flow such as K = 0.25 
and obtain a fast settling but low, level at the destination; or else liberally in- 
crease K to the allowed maximum value of 0.35 for stability, causing level 
fluctuations which take a longer time to settle but obtaining a higher final 
level. In summary, the graph in Figure 18 reveals an unstable system, for 
K = 0.45. For K = 0.25, one has a definite stable case. And for K = 0.35, the 
system, though stable, is very close to be uncontrollable. The simulation is 
therefore in complete agreement with the predictions of the Satche technique 
applied in the two-delay case. 


D. OTHER NON- TRIVIAL PURE TIME DELAY CASES: 

Since the computer program is a general one, numerous delay cases can 
be examined. 

(i) Single Pure Time Delay: 

This situation is the simplified version with no time lag illustrated in 
Appendix D, (Equation D-5). Figures 21 through 23 present three cases for 
K = 2, 1 , 77/2 and which are unstable and optimal, respectively. In the last case 
77/2 has been approximated by 1.57. The Satche functions F, (a;) and F 2 (tj) for 
equation (D-5) are i", a straight line through the imaginary axis and -Ke~ ia ’, a 
circle of radius K. Furthermore, the fact that F 2 in Figures 21 and 23 looks 
like an ellipse rather than a circle is no surprise since it is simply due to the 
graphic subroutine trying to optimize the area where the diagram is to be drawn; 
numerically, Figures 21 and 23 represent a circle. Figures 21 through 23 are 
in agreement with the results in Refs. 42 and 43. 

Figures 24 and 25 show an unstable and a stable single pure time delay 
case with time lag, respectively. The characteristic equation for these two 
cases is taken as equation (57) with b 2 = 0. These two figures illustrate, 
somewhat, how to stabilize a process by manipulating K and v . A decrease in 
7 " l in Figure 25 flattens to the right the parabola of Fig. 24 as dictated by 
equations (64-65). This allows a larger region of Fi to the left of F 2 to be ex- 
posed and on which the critical part of F 2 '(***) can possibly fall, thus producing 
a stable case. Another manipulation involves a decrease in K thus concentrat- 
ing the critical part of F 2 (***) to the left of F } , as seen in going from Figure 
24 to 25. Here, one also notes that a decrease in the time delay improves the 
stability. 
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(ii) Extension to Cases with Three, Four and Five non-trivial 
Pure Time Delays: 

The preceding model is extended to incorporate up to five pure time 
delays. The general equation (43a) and Figure 4d are utilized with N = 3, 4, and 5. 
The characteristic equation for three delays is an extension of equation (57). It 
is given as: 

F(s, e s ) = s(T lS + 1) + bjKe- 81 ! + b 2 Ke" s(t i + l 2> + b 3 Ke~ s( tl+t2+t3) = 0 (75) 

with 

Fj(s) = s(r 1 s + 1) (76) 

and 

F 2 (e s ) = - [b 1 Ke -st 1 + b 2 Ke" s( 1 i +t 2) + b 3 Ke- s ( t i +t 2+ t 3)] (77) 

For four delays equation (57) becomes: 

F(s , e s ) = s(r 1 s + 1) + bjtfe -stl 

+ b 2 Ke“ , <‘i +t a ) + b 3 Ke- s < t i +t 2 +t 3> (78) 

+ b 4 Ke- s(t i +t 2+ t 3+ t 4) = 0 

with 

Fj(s) = s(Tj s + 1) (79) 

and 

F 2 (e s ) = _ [bjKe- 81 ! + b 2 Ke- s ^i +t 2> (80) 

+ b 3 Ke- s ( t i+t2+t3+t 4 )] 
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And for five delays, one has: 


F(s, 


e s ) = s(TjS + 1) + bjKe -511 


+ b 2 Ke -s ( tl+t2 ^ 


+ b 3 Ke- s(t i +t 2 +t 3) + b 4 Ke- s(t i +t 2 +t 3+*4) (81) 

+ b s Ke- s(t l +t 2 +t 3 +t 4 +t 5)] 


with 


Fj(s) = s(t 1 s + 1) 


(82) 


and 


F 2 (e s ) = - [b^e- 811 + b 2 Ke- s(tl + t2 > + b 3 Ke " s < t l+ t 2 +t 3 ) 


b Ke -s(tl + t2+t 3 +t 4) +b Ke _s<tl + t 2 +t 3 +t 4+ t 5>] 


(83) 


Figures 26 through 28 illustrate near optimum cases (optimal K) for the three - 
pure time delay problem. Figures 29, 30, 31 and 32 presents diagrams for the 
case of four and five pure time delays, respectively. Figures 29 and 30 shows a 
stable situation while Figure 31 indicates instability. Figure 32 gives a stable 
case for five pure time delays. 
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SECTION VI 
CONCLUSION 


A. JUSTIFICATION OF THE DIFFERENTIAL-DIFFERENCE EQUATION: 

The models forged so fax, containing a combination of continuous and dis- 
crete equations (with derivatives and differences) can be simplified to a purely 
discrete form involving only difference equations. Such a discrete model poses 
no difficulties in the stability analysis. In this case, the characteristic poly- 
nomials involve the Z transform which acts efficiently on transcendental terms. 
As seen, this is not so in the case of the continuous-discrete formulation. In 
justification of the continuous approach taken here, the adjectives discrete and 
continuous are relative terms depending on the application. As an illustration, 
feedback information flow entering the decision-making headquarters in the 
form of a telephone call once a day, for instance, appears as a discrete event 
among other daily office activities. However, if these daily phone reports serve 
to inform the management on a transportation network suffering time delays of 
monthly magnitude, then they form a continuous stream while the transportation 
flow is of much more discrete nature. 


B. COMPUTER TIME: 

The total computer time for analyzing the stability of the models presented 
so far requires at most 30 seconds of processing per case. Therefore the con- 
templation of developing more efficient routines is not pursued. 


C. THE CHOICE BETWEEN THE LAPLACE APPROACH AND THE 
TIME DOMAIN APPROACH: 

The Laplace transform method used here differs noticeably from the theory 
developed in the time domain, often called the state variable theory. The first 
one is rigorous since it optimizes a model by insuring its stability. The second 
approach, instead, heuristically optimizes an objective function. This function 
is often left at the analyst discretion to formulate, based on what he feels are 
important variables. 


Preceding page blank 
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D. SUMMARY 


, R. Bellman 53 in a survey on the literature related to Kalecki's Model 22 , 
discusses two characteristic equations, one of the form: 


A.s 2 - s + b - ce -s = 0 


( 84 ) 


where \ b and c are constants, and the other one as: 


s = b(e _stl - e _st2 ) 


(85) 


where b is a constant and t 2 = t 1 +1. The above author comments that the 
determination of the roots of eq. (84), the one containing a simple delay, has 
been fully examined for stability. Equation (85) on the other hand, containing 
two coupled delays is, quoting the author, "considerably more difficult” and 
where 'the details are much more perplexing" in the analytical stability ap- 
proach. The extended Satche method developed herein, can easily determine 
the stability of eq. (85) from the graphical point of view, thus by-passing the 
analytical determination of the roots. In addition, since this technique is 
computer-based, it serves as a fast and accurate method for a particular user. 

Furthermore, the extension to multiple delays is quite impossible without a 
computer, justifying the use of a machine. The complexity of Figure 31 supports 
such statement. 

One notes that, in general, as the number of pure time delays increases, the 
allowable value for K decreases allowing for less resource flow. The undesir- 
ability of delays in the real world is thus mathematically illustrated. Most im- 
portant is the total agreement between the simulation and the criterion predic- 
tion. This Satche extension must undergo other tests if it is to become an es- 
tablished technique. Finally, the validity of the model for multiple inputs and 
outputs and for non-deterministic random variables remains to be demonstrated, 
thus providing grounds for further research. 


34 



BIBLIOGRAPHY 


1. Bernoulli, J., "Meditationes. Dechordis Vibrantis ...» Commentari 
Academia Scientiarum Imperialis" Petropolitanae, 3 (1728), 13-28, Collected 
Works, v. iii, p. 198. 

2. Lagrange, J. L., "Sur 1* integration d* une equation differentielle a differ- 
ences finies, qui contient des suites recurrentes," Miscellanea Philosophic- 
Mathematica Societatis Privatae Taurinensis, 1, (1759). 

3. Poisson, "Sur les equations aux differences melees," J. Ecole Polytech., 

Vol. 6, cahier 13, (1806), pp. 126-147. 

4. Cauchy, A., "Memoire sur l'emploi des equations symboliques dans le 
calcul infinitesimal et dans le calcul aux differences finies," C. R. Acad. 

Sci. Paris., Vol. 17, (1843), pp. 449-458. 

5. Bellman, R., Cooke, K. L., "Differential- Difference Equations," Academic 
Press, (1963). 

6. Pinney, E., "Ordinary Difference- Differential Equations," U. of Calif. Press, 
(1959). 

7. Ogutzoreli, M. N., "Time-Lag Control Systems," Academic Press, (1966). 

8. Wright, E. M., "Linear difference-differential equations," Quart. J. of Math., 
Vol. 17, p. 245, (1946). 

9. Bateman, H., "The Solution of a System of Differential Equations occurring 
in the Theory of Radio-active Transformations," Proc. of the Cambridge 
Philos. Soc., Vol. 15, (1910), pp. 423-427. 

10. Gerasimov, S. G., "Fizicheskie osnovy dinamiki regulirovaniia teplovykh 
protsessov" (The Physical Bases of the Dynamics of Regulation of Heat 
Processes), Energoizd, Moscow, (1937), pp. 71-75. 

11. Lagrange, J. L., "Recherches sur la nature et la propagation du son," 
Miscellanea Philosophico-Mathematica Societatis Privatae Taurinensis, 1, 
(1759). 


35 



12. Arley, N., "On the Theory of Stochastic Processes and Their Application 
to the Theory of Cosmic Radiation," G.E.C. Gads Forlag, Copenhagen, 

(1943). 

13. Minorsky, N., "On Certain Applications of Difference- Differential Equations," 
Dept, of Engineering and Mathematics, Stanford U., (1948), 38 pages. 

14. Nisolle, L., "Sur la stability des regulations a impulsions retardees ou 
amorties," C. R. Acad. Sci. Paris., Vol. 211, (1940), pp. 762-765. 

15. Volterra, E., "On Elastic Continua with Hereditary Characteristics," J. 
Applied Mech., Vol. 18, (1951), pp. 273-279. 

16. Picard, E., "La mecanique classique et approximations successives," Riv. 
di Sci., Vol. 1, (1907). 

17. Satche, M., J. Applied Mechanics, Vol. 16, p. 419, (1949). 

18. Bellman, R., Harris, T. E., "On the Theory of Age-dependent Stochastic 
Branching Processes," Proc. Nat. Acad. Sci., Vol. 34, (1948), pp. 601-604. 

19. Lotka, A. J., "Relation between Birth Rates and Death Rates," Science, 

Vol. 26, (1907), p. 21. 

20. Haldane, J. B. S., "A Mathematical Theory of Natural and Artificial Selec- 
tion," Proc. Cambridge Philos. Soc., Vol. 23, (1927), pp. 607-631. 

21. Volterra, V., "Lemons sur la theorie mathematique de la lutte pour la vie," 
Gauthier-Villars & Cie., Paris, (1931). 

22. Frish, R., Holme, H., "The Characteristic Solutions of a Mixed Difference 
and Differential Equation occurring in Economic Dynamics," Econometrica, 

3, (1935), pp. 225-239. 

23. Kalecki, M., "A Macrodynamic Theory of Business Cycles," Econometrica, 
Vol. 3, (1935), pp. 327-344. 

24. James, R. W., Belz, M. H., "The Significance of the Characteristic Solutions 
of Mixed Difference and Differential Equations," Econometrica, Vol. 6, 

(1938), pp. 326-343. 

25. Samuels on, P., "Foundation of Economic Analysis," Harvard U. Press, 
Cambridge, Mass. (1947). 


36 



26. Bateman, H., "An Integral Equation Occurring in a Mathematical Theory 
of Retail Trade," Messenger of Math., Vol'. 49, (1920), 134-137. 

27. Theiss, E., "Dynamics of Savings and Investment," Econometrica, Vol. 3, 
(1935), pp. 213-224. 

28. Tinbergen, J., "Annual Survey: Suggestions on Quantitative Business Cycle 
Theory," Econometrica, Vol. 3, (1935), pp. 241-308. 

29. Tsien, H. S., "Engineering Cybernetics," McGraw-Hill Book Co., N.Y., 
(1954). 

30. Smith, O. J. M., "Feedback Control Systems," Chap. 10, (1958) N.Y.: 
McGraw-Hill. 

31. Truxal, J. G., "Feedback Theory and Control System Synthesis," N.Y.: 
McGraw-Hill, (1955). 

32. Rogers and Connolly, "Analog Computation in Engineering Design," McGraw- 
Hill, (1960). 

33. Flagle, C. D., Huggins, W. H., "System Dynamics," "Operations Research 
and Systems Engineering," John Hopkins Press, Baltimore (1960). 

34. Oetker, R., "On the Control of Sectors with Dead Time," Proceedings of 
First International Congress of the Int. Fed. of Auto. Control, Moscow, 

1960, Vol. 1, Butterworth (1961). 

35. Paxil, J., "Analog Simulation," San Fernando State College, Calif., Course 
taught at the Civil Service Commission, Washington, D.C. (1972). 

36. Coughanour, D., andKoppel, L., "Process Systems Analysis and Control," 
N.Y.: McGraw-Hill, (1965). 

37. Tyner, M., May, F., "Process Engineering Control," N.Y.: Ronald Press, 
(1968). 

38. McNaughton, R., and Papert, S., "Counter- Free Automata," MIT Press, 
Cambridge, Mass., (1971), pp. 61-66. 

39. Forrester, J. W., "Industrial Dynamics," MIT Press and J. Wiley and Sons, 
N.Y., (1961). 


37 



40. Gue, R. L., Thomas, M. E., "Mathematical Methods in Operations Research," 
The MacMillan Co., N.Y., (1968). 

41. Prabhu, N. V., "Queues and Inventories, A Study of Their Basic Theory," 
Wiley, N.Y., (1965). 

42. Stoker, B., "Stability Criteria for Linear Dynamical Systems," N.Y.: 
Academic Press, (1968). 

43. Distefano, Stubberud, and Williams, "Feedback and Control Systems," 
McGraw-Hill, (1967). 

44. Choksy, N. H., "Time Lag Systems — A Bibliography," IRE Trans. Autom. 
Contr., Vol. AC-5, No. 1, Jan (1960), pp. 66-70. 

45. Single, C. H., and Stubbs, G. S., "Transport Delay Simulation Circuits," 
Westinghouse Electric Corp., Atomic Power Division, Rept. WAPD-T-38, 
(1956). 

46. Krall, A. M., "An Extension and Proof of the Root Locus Method," J. 

SIAM 9 (1961), pp. 644-653. 

47. , "A Closed Expression for the Root- Locus Method," J. SIAM, 

11, (1963), pp. 700-704. 

48. , "An Algorithm for Generating Root-Locus Diagrams," Com- 

munic. of the ACM, Vol. 10, No. 3, March (1967), pp. 186-188. 

49. Chu, Y., "Feedback Control System with Dead- Time Lag or Distributed Lag 
by Root- Locus Method," Trans. Amer. Inst. Elect. Engrs., Vol. 71, 291, 
(1952), pp. 291-296. 

50. Yuan-Xun, Q., Iong-Quing, L., and Li an, U., "Effect of Time-Lags on Sta- 
bility of Dynamical Systems," Proceedings of First International Congress 
of the Int. Fed. of Autom. Control, Moscow, (1960), Vol. 1, Butterworth, 
(1961). 

51. Shaughnessy, J. D., and Kashiwagi, Y., "The Determination of a Stability 
Indicative Function for Linear Systems with Multiple Delays," NASA Tech. 
Rpt. TR R-301, (1969). 

52. IBM, "System/360 Continuous System Modeling Program," GM20-0367-2, 
(1969). 


38 



53. Bellman, R., "A Survey of the Mathematical Theory of Time-Lag, Retarded 
Control, and Hereditary Processes,” Project Rand, R-256, (1954), The Rand 
Corp., Santa Monica, Calif., p. 88. 


39 



Receding page blank: not femet? 


APPENDIX A 

BASIC DEFINITIONS USED IN THIS WORK 


1. A system or process: an arrangement of elements connected or related in 
such a m ann er as to form a working entirety. 

2. A control system: a system able to regulate or command itself or other 
systems. 

3. A feedback control system: a control system where the output of some 
controlled system variable is compared with the input to the system so that 
the appropriate control action may be formed as a function of the output or 
input. Feedback increases the accuracy of a control system and favors 
overall system stability. 

4. Modeling: a mathematical method used to express the real world through: 

(a) Functional or differential equations. 

(b) Block diagrams, similar to flow charting. 

(c) Flow graphs. 

Block diagrams and flow graphs are schematic representations of the model's 
equations and real components. Appendix B explains the block diagram notation. 
The flow graph method appears in Appendix C. 
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APPENDIX B 

TERMINOLOGY OF FEEDBACK AND CONTROL: THE BLOCK 
DIAGRAM AND THE TRANSFER FUNCTION 


A. THE BLOCK DIAGRAM: 

The block diagram notation is a graphical representation of the overall 
cause and effect relationship between input and output of a process or system. 
The components or elements of a system is each characterized by a box (or 
block) with an arrow penetrating and exiting, representing an input and output, 
respectively. The box usually contains a mathematical transformation such as a 
differentiation, integration, division or multiplication which acts on the input to 
produce the desired output. Figure la illustrates such a representation. All 
the blocks are related to one another and eventually form the overall system. 
Figure lb illustrates a generalized beedback control system. In this diagram, 
the control block function can be varied, the process is fixed and the feedback 
is subtracted (negative feedback) or added (positive feedback) from an input 
originating outside the system. 


B. THE TRANSFER FUNCTION (T.F.): 

The ratio of the output to the input of a system constitutes the transfer 
function. This ratio can be determined analytically by equating the output to the 
transformed input for each block of the system. This gives a series of coupled 
equations which can be combined into a single equation containing the original 
input and final output. The ratio of these two is then readily obtained. 

The method in Appendix C provides a more elegant and rigorous way of 
determining the transfer function for linear, deterministic systems. 
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APPENDIX C 

NETWORK FLOW GRAPHS 
AND MASON'S LOOP RULE 

A. NETWORK OR SIGNAL FLOW GRAPHS: 

Basically, the block diagram notation in Figure lb can be simplified to an 
algebraic graph as follows: 


INPUT 


H=NET FORWARD FLOW 
_G=NET FEEDBACK FLOW 



where the transfer function is 


OUTPUT _ T F _ GH 
INPUT ‘ " 1 - GH' 


The arrows in the above diagram are called branches while the destination dots 
are called nodes. A path is a branch or a continuous sequence of branches which 
can be traversed from one signal (node) to another signal (node). A loop is a 
closed path which originates and terminates on the same node and,' along which, 
no node is met twice. 


B. MASON’S LOOP RULE: 

Equation (C-l) is obtained by equating the incoming flow of values to the 
outgoing one at each node, since no accumulation can occur in a flow graph. 
This produces a set of equations which reduce to equation (C-l). A shorter 
way to derive the transfer function is to utilize Mason's Loop Rule (see ref. 
43). The Rule states that, given a flow graph, the transfer function between 
two nodes, i and j respectively, is given by: 
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( 02 ) 


T. 


F. 


i j 



A. . 


i j 


k 


A 


where 

Pj . = k th forward path from variable nodes i to j, 

A = determinant of the graph, 

A. . = cofactor of the path P. . 

^k 

and the summation is taken over all possible k paths from the k to j nodes. The 
cofactor Ai j k is the determinant with the loops touching the k th path removed. 
The determinant A is: 


M, Q 


A = 1 


-I>» + L L - L - -H 


L L . L _ + 


(03) 


m= 1 , q= 1 


where L q equals the value of the q th loop transmittance. Therefore, in words, 
equation (C-3) signifies that A = l - (sum of all different loop branches) + (sum 
of branch products of all combinations of two non-touching loops) - (sum of the 
branch products of all combinations of three non-touching loops) + . . . . Two 
loops are non- touching if they do not have a common node. 


C. ILLUSTRATION OF MASON'S LOOP RULE: 
Consider the following flow graph: 
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There are five forward paths connecting the input to the output: 


p ! =H, 

h 2 h 3 h 4 

(03) 

p 2 = Hi 

g i H 3 G 3 

(04) 

p 3 = H 1 

G 1 H 3 H 4 

(05) 

P 4 = H 1 

h 2 h 3 g 3 

(06) 

P 5 

~ P 5 

(07) 


There are four self loops: 


Li = H i H 2 G 2 

(C-8) 

L 2 _ P 4 G 4 

(09) 

L 3 = H x H 2 H 3 G s 

(O10) 

L 4 = H 2 H 3 G 6 

(Oil) 


Loops Lj and L 2 do not touch. Therefore, the determinant is: 

A = 1 - (Lj + L 2 + L 3 + L 4 ) + (L t L 2 ) (C-12) 

The cofactor along P x is obtained by removing from A the loops that touch P : 


Aj = 1 - 0 = 1 


(013) 


* 
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The same applies for P 2 , P 3 , P 4 and P s : 


A 2 = 1 - 0 = 1 

(014) 

A 3 = 1 - 0 = 1 

(C-15) 

i-H 

II 

O 

1 

rH 

II 

<3 

(C-16) 

*5 = 1 - L « 

(C-17) 

and therefore, 


j p _ Aj + P 2 A 2 + P 3 A 3 + P 4 A 4 + P 5 A s 



A 

- (H, H 2 H 3 H 4 ) + Hl G x H 3 G 3 + H| G x H 3 H 4 + H g (1 - H, H 3 G fi ) 

1 - (H, H 2 G 2 + H 4 G 4 + H x H 2 H 3 G s + H 2 H 3 G 6 ) + H, H 2 H, G 2 G, 

(C— 18) 
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APPENDIX D 


THE NYQUIST CRITERION 

A. CONFORMAL MAPPING, THE ENCIRCLEMENT THEOREM: 

If the variable s, in the Laplace domain, assumes real and imaginary parts 
as follows: 

s = cr + i co (D-l) 


(where i = /-l), then the characteristic equation in s also assumes real and 
imaginary parts given by: 


F (s, e s ) = U (cr, co) + i V (cr, co) 


(D-2) 


According to the theory of complex variables, a conformal mapping from s to 
F is defined as follows. To any point (cr, co) in the s-plane, there corresponds 
a unique point (U, V) in the F-plane. If the point (cr , co) traverses a closed con- 
tour C in s then (U, V) will trace a corresponding contour, F . Figure 6a illus- 
trates such mapping. Accordingly, the Encirclement Theorem states that when 
C is traversed once clockwise by a point (cr, co) then the corresponding mapped 
point (U, V) revolves N times counter-clockwise on F in Figure 6a, such that, 


N = - Z 


(D-3) 


where Z is the number of those roots of equation (D-2) which lie within C pro- 
vided that none of them lie on C itself. A negative N indicates that F is tra- 
versed in clockwise fashion, while a positive N indicates a counter-clockwise 
direction. A complete proof is available in Reference 42. 
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B. THE CONDITION ON THE s-PLANE FOR ASYMPTOTIC STABILITY 


The condition for asymptotic stability, exposed in Section IVA, requires that 
all the roots of the characteristic equation to have negative real parts and there- 
fore to lie in the left-half of the s-plane in Figure 6a. This is briefly explained 
as follows. The response of a linear stable system to an impulse function must 
tend to zero as t - as stated in Section IVA. This response is generalized as: 


y (t) = Aj e 


(ai + iwj ) t 


+ A„ 


(0- 2 + i£ l ) 2 ) t 


( cr n + i£ ‘’n) t 


(D-4) 


(i = v^T) 


where A^ 9 A ^ • • • A^ sire coustduts dud (o - ^ 4" 1 co^ )» (o- 2 + i ^ 2 )> • • • * ( a n + ia, n ) 
are the complex roots of the characteristic equation (D-2). It can be seen from 
equation (D-3) that if the cr ’ s are negative, then y(t) decays exponentially thus 
satisfying the requirement for asymptotic stability. Abiding to such prerequisite, 
one cancels all possibilities for unstable roots to lie in the left-hand plane. This 
leaves only the co -axis and the right-hand s-plane in Figure 6a to examine for 
unstable roots. The general contour C then specializes to a newcontour, C', 
which is shown in the left portion of Figure 6b. It consists of the whole imaginary 
axis between co =~ cc to w = + co , along with the infinite semi-circular arc en- 
compassing an area to the right of the direction of the arrow, namely, the com- 
plete, bounded right-hand plane. 

The conformal mapping of this semi-circular closed contour C' onto the 
F-plane produces a closed T contour. First, one considers the semi-circular 
arc portion (excluding the co axis) of C* on which s can be written in terms of 
polar coordinates as: 


s = Re i9 (tt/ 2 >6 > - tt/ 2) (D-5) 

where R = infinity, i = /-l. The mapping of this portion onto the F-plane gene- 
rates a clockwise V portion which also bounds the entire right hand F-plane 
with the imaginary F axis. The equation of this r portion is immaterial since 
every point on it is at an infinite distance from the origin. It can be taken as a 
semi-circle with infinite radius. It is referred to as part of the F-locus at in- 
finity. Returning to the s-plane, the imaginary s-axis remains to be examined 
from s = gj = -co to s = w = +co. Depending on F, the mapping produces various 
loci which then connect with the part of the F-locus at infinity to form a closed 
contour. The area enclosed depends again on the direction of the arrows on r. 


50 



C. THE NYQUIST CRITERION: 

The Criterion basically combines the requirements of section (A) with those 
of section (B) by determining the number of roots which reside inside C' of the 
s-plane in Figure 6b, (i.e. those roots with positive real parts and therefore 
those which cause instability). It can be stated as follows: 

According to the Encirclement Theorem, the number of times, N, that the F- 
locus encircles the origin in a clockwise sense corresponds to the number of 
roots within C' in the s-plane, i.e. those roots with positive real parts and which 
cause the system to be an unstable one, according to the stability condition on 
the s-plane. Therefore, if the F-locus never encircles the origin completely, 
then the system is (asymptotically) stable, and unstable otherwise. 


D. ILLUSTRATION: 

The previous criterion is best illustrated as follows. Consider the charac- 
teristic equation containing one pure time delay as : 

F (s, e s ) = s + K e“ s = 0 (D-6) 


where K is an arbitrary constant. 

The mapping of s onto F(s, e s ) consists of two parts. First, the semi- 
circular portion of C* in the s-plane is represented by Equation (D-4). As 
previously stated, the F(s, e s ) contour corresponding to every point (o-, co) on 
that section is itself a semi-circular curve traced at an infinite distance from 
the origin. It is shown in the F-plane drawn in the right portion of Figure 6b 
and labeled | F(s, e s ) | - 1X5 . The second part of the mapping traces F (iw, e 1<w ) 
for w varying from - “to +oo. This tracing is done by firstly” substituting ico for 
s in equation (D-5). The result is 


i <*> + K e l “ = 0 


(D-7) 


Since e 1 " = cos co - i sin co , one separates the real and imaginary parts of 
equation (D-6) to yield: 
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Xj (co) = Real F = K. [cos w] 


(D-8) 


Yj (co) = Imaginary F = K [i (co - sin co)] 


(D-9) 


For every won the imaginary axis of the s-plane, there corresponds an ordered 
pair (Xj (w), Yj (co)) on the F-plane. Numerical values are calculated for X l 
and Yj for a varying co . The result is the mapping of co onto F. It is shown as 
the "wiggle-like" portion along the imaginary axis of the F-contour in the right- 
hand portion of Figure 6b. One notes that three curves appear each correspond- 
ing to a different value of K, namely 1.0, 77/2, and 2.0, respectively. All three 
curves connect with the semi-circular portion of F labeled |F(s, e s )| = 00 at 
infinity. The mapping, called the Nyquist plot, therefore completely encloses 
an area to the right of the direction of the arrows (shaded in Figure 6b) in F. 

It can be seen that for K = 2 the F-locus encloses the origin twice in the 
clockwise sense, once very close to it and another time very far from it at 
R = o° so that N = -2 and therefore Z = -N = 2, indicating two roots with positive 
real parts in C*. The cases for K = 1.0 and 77/2 on the other hand reveal no 
encirclement and therefore constitute stable cases. The value K = 77/2 is at the 
threshold of stability. It represents the largest allowed value for stability. 
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APPENDIX E 


LISTING OF THE COMPUTER PROGRAM DEVELOPED TO IMPLEMENT 
THE SATCHE CRITERION FOR MULTIPLE DELAYS AND TO PERFORM 

COMPUTER GRAPHICS 


A. THE SUBROUTINES USED: 

The following program is written in FORTRAN IV . The implementation 
utilized an IBM 360/91 operating under Release 20, a relatively high speed, 
large memory machine. 

Subroutine SCAN determines the intersection points while subroutine 
SATCHE converts the numerical results into a special code which is transferred 
on a 9-track magnetic tape. This tape is then mounted on a CALCOMP 780 
(California Computing) System which in turn draws the desired plots on a drum 
using ink and paper. 

The subroutines PLOTS, PLOT, LINE, NUMBER, SCALE, and AXIS are 
made available by the CALCOMP company and assist in drawing the plots ac- 
cording to a desired scale, format and size. CALCOMP also provides a 
repertory of symbols, letters, and other characters which label the 
diagram. 

B. LIST OF COMPUTER SYMBOLS IN CONTEXT WITH THE MAIN 

WORK (IN APPROXIMATE ORDER OF APPEARANCE): 


XI = (co) = Re Fj 
X2 = X 2 (co) = Re F 2 

Y1 = Yj (co) = Im Fj 
Y2 = Y 2 (co) = Im F 2 

BUFFER = A dummy vector used as a "scratch pad" by subroutine 
SATCHE. 


53 



OMEGA. 

B 

T 

PI 

PIB 

DPIB 

I 

NFMAX 

NFILE 

M 

XK 

TAU1 

MM 

Nl, N2, N3, 
W 

XDPI 

XT 

XX2 

YY2 


= a vector containing all the values for co. 

= a vector containing the feedback factors, b i} in order of ap- 
pearance in the network forward flow. 

= a vector containing the exponents of the delays in order of ap- 
pearance in the network forward flow. 

= 7 T 

= selected maximum value for co. 

= selected increment value for co . 

= subscript for OMEGA, XI Yl, X2, Y2 initially set at 1 and 
finally set at 513. 

= selected maximum number of computer runs per batch. 

(one run per delay case) 

= file number assigned to each case, initially set at 1 and finally 
set at NFMAX. 

= number of pure time delays in the particular problem being 
processed. 

= K = management control 
= t = time lag 
= dummy subscript 

N1MAX, N2MAX, N3MAX = various pointers of secondary im- 
portance used to produce a convenient output format. 

= an intermediary value for OMEGA. 

= an increment smaller than DPIB (of secondary importance) 

= a temporary value for T 
= a temporary value for X2 initially set at 0 
= a temporary value for Y2 initially set at 0 
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K 

TEMPX2 

TEMPY2 

WP 

N 

IMAX 

JMAX 


= a dummy subscript 
= another temporary value for X2 
= another temporary value for Y2 

= a value used to transform co into a fractional part of 77 to be- 
come readable in the output printout. 

= a code indicating "intersection" or "no intersection." 

= 513 = maximum number of values for XI, Y1 and of OMEGA 
used to evaluate XI and Yl. 

= 513 = maximum number of values for X2, Y2 and of OMEGA 
used to evaluate X2 and Y2. 


IIMIN = subscript of the OMEGA corresponding to co y , the smaller 
value at which one of the outermost intersections occurs. 


UMAX = subscript of the OMEGA corresponding to <^ 2 the larger value 
at which the other outermost intersection occurs. 

J = temporary subscript 

ID = the number of values of OMEGA for Fl, XI and Yl on FI 

between the outermost intersections. 

IIIMIN = subscript of a reconstructed OMEGA for Fl, XI and Yl with 
the subscript ID = 1 to ID = UMAX - IIMIN. This is necessary 
so that the curve Fl and F2 have the same order of magnitude 
when the Satche diagram is performed. 


SUBROUTINE SCAN: (additional symbols to the above) 

D = distance between (XI, Yl) and (X2, Y2). 

DMIN = smallest distance between (XI, Yl) and (X2, Y2). 

JX = subscript identifying the point (X2, Y2) closest to the given 

point (XI, Yl). 
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SUBROUTINE SATCHE: (additional symbols to the above) 


XL, YL 
XSUBX, 

D 

JJ, NN 
XN 


= values specifying the paper area on which the Satche diagram 
is to be drawn. 

YSUBY, YSUBX, XSUBY = CALCULATED origin of the Satche dia- 
gram. (Optimized values such that the diagram is centered 
within the area specified by XL, YL). 

= a temporary value used in the above optimization. 

= dummy subscripts. 

= a value used for proper format (of secondary importance) 
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nno o n onnnn on 


VN IV G LEVEL 20 


MAIN 


DATE 


7221 4 


16 / 08/18 


C 

C NYQUIST-SA TCHE DIAGRAM 

DIMENSION XI! 600) , X2 ( 600) , Y 1 ( 60 0 )» Y2 ( 60 0 ) . 3UFF = P ( 750 0 ), OMEGA ( 60 0) • 
1 B ( 10) tT( 10) 

C 

c 

C SETTING UP A WORK AREA FOR PLOTTING 

CALL PLOTS! BUFFER * 7500,0) 

C 

C 

c 

c 

C CONSTANTS DEFINED 

P 1 = 3, 141 59 
P IB— P I 
DP 1=16.416. 


reading the number of computer runs DES I R £ ) 

READ (5.77) NFMAX 
77 FORMAT (12) 


DO 888 NFILE=1 , NFMAX 

READING IN ORDER THE DESIRED NUMBER OF RUNS, 

FEEDBACK MULTIPLICATIVE FACTORS, AND THE EXPONENTS OF THE 
PURE TIME DELAYS. 

READ ( 5,99)NFILE.M ,XK,TAU1 ,(B(MM),MM=1 ,5> ,(T(M*> » MM= 1 ,5) 

99 FORMAT ( 12. I1*12(F5.2>) 

1=1 

771 WRITE! 6, 1 )NF ILE.M .XK.TAUl , < 9 ( M M > , MM= 1 . 5 ) » (T(MV ) • MM=1 .5 ) 

1 FORMAT ( • 1 • • /// • 3X , 

2 1 COM PUTER RUN NUMBER *,I2. 9X,*SATCHE D A G EN cR AT I ON • / • 3X , 
7 • THE NUMBER OF DELAYS IS *,I2./.3X. 

8* THE CONSTANTS ARE: # ,/3X, 

1 • K= 1 ,F8.3, 5X, • TA U1 = * ,F8.2, // .3X. 

3' B( 1 ) TO 0 ( M) : • »5( F8 .3,2 X) , /. 3X • 

S'TIl) TO T ( M ) J • ,5(F8.3,2X)> 

772 WR I TE (6,3) 

3 FORMAT!//, 4X, ‘OMEGA OF F1*,8X.*FE F1S10X, 

1*IM FI* , 1 0 X . *RE F2*,10X* # IM F2‘, 9X,'CMEGA OF F2*> 

THE NEXT TWO STATE ME NT S ARE USED FOR A 
DESIRED PR INT-FORM AT . 

N 1= 1 

NlMAX=26-< M*2) 

773 W=-PI 
XDP 1=! P I/DP I )/2. 

C GENERATION OF F(S). 

4 XI! I >=-TAUl*W**2 

61 Y1(I) = W 

C GENERATION OF F ( E XP( S ) ) , THE CURVE WITH DELAYS. 

XT=T< 1 ) 

X X2= 0 • 
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MAIN 


DATE = 7221 * 16/08/18 


YY2=0 . 

DO 67 K = 1 , M 

TEMPX2=B( K )*XK*COS (XT*«) 

TEMP Y 2=8 ( K)*XK*SIN(XT*W) 

IF(K-M) 69 ,68.68 
69 XT=XT + T ( K + 1 ) 

68 X2! I >=- ( XX2+TEMPX2 > 

XX2 =-X2( I ) 

Y2( I )= Y Y 2+ TEMPY 2 
YY2=Y2< l ) 

67 CONTINUE 

WP=256.*W/PI 
OMEGA! I )=WP 

777 WRITE(6.2) WP • X l ( I ) • Y1 ( I ) • X2 ( I ) • Y2 C I ) . WP 

2 FORMAT ( 3X » F 6* 1 , */2 56 PI *, 5 X .4 ( 1 P El 0 . 3 . 5 X) * OPF 6 . 1 , */ 256 PI*) 

IF(N1— N1MAX)8. 10* 1 C 
10 N 1 = 0 

WRITE! 6 . X ) NF I LE , M .XK.TAUl * ( 8 (MM ) , MM= 1 *5) , (T ( MM) , MM= 1 ,5 ) 

WRITE! 6.3) 

8 N1=NU1 

88 I F( ( ,W+ XDP I )— P 1 8 ) 1 2 *9 .9 
12 1=1+1 
7 w = W+P I/DPI 
GO TO 4 

9 CONTINUE 

19 WRITEC6.1 )NFILE .M .XK.TAUl • (B (MM ),MM = 1 ,5) ,(T ( NN) . MM=1 .5 ) 

97 WR I TE ( 6 * 23 ) 

23 FORM AT ( / * 3 X . • F I ND t NG ALL POINTS CF INTERSECTION') 

96 WR I TE { 6 * 20") 

20 FORMAT! /. 3X . *OBSER VE THE FOLLOWING CODES BELOW • ,/ 3X • 

1 * N= 0 2 NOT A REGION OF I NTE RSEC T 10 N * . / , 3X . 

2 * N= 1 : A REGION OF INTERSECTION*) 

95 WRITE! 6.21 ) 

21 FORMAT! //,4X, * OME GA OF F1».6X.*PE F 1 * * 8X . 

1*IM FIS Q X • 'RE F2*.08X.*IM F2» . 7X. ‘OMEGA 3F F2*.4X.*N*> 

98 IMAX=I 
JMAX=I 

I IMA X=0 
I IMIN=999 
N 2= l 

N 2M A X= 1 9— ( M*2 ) 

DO 32 1=1 , IMAX 

CALL SCAN! I.Xl(I).Yl(I).X2. Y2 *N . L. JMAX) 

24 WR I TE ( 6 . 26 ) OMEGA! l).Xl(I).Ylfl) *X2 ( J) « Y2 ( J) , C MEG M J ) .N 

26 FORMAT! 3X.F6. 1 . */2 56 PI * . 3 X • 4 < l P El 0 . 3 .3 X) . ) PF 6 . 1 • •/ 256 PI*.2X.I2) 
IF(N2-N2MAX)30.31 , 31 

31 N2=0 

WRITE! 6. 1 )NFILE.M .XK.TAUl . (B(MM ).MM = 1 .5) . (T!*M.) • MM=1 .5) 

WRITE (6.23) 

WRITE! 6,20 ) 

WRITE (6.21) 

30 N 2=N2f 1 

C LOCATING THE OUTERMOST INTERSECTIONS. 

45 IF(N) 3 2.32.27- 

27 IF( I- I IM IN >29.2 9.3 4 
29 I IM IN= I 
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MAI N 


DATE 


72214 


16/08/ 18 


J 1 = J 

34 IF(I-IIMAX) 32*33,33 
33 IIMAX=I 

J2=J 

32 CONTINUE 

WR ITE( 6* l ) NF I LE *M ,XK,TAUl * C 9 ( M M ) » MM = 1 .5) . <T ( MM) » MM=1 *5 ) 
WRITE(6,35> 

35 FORMAT! //. 3X, • THE OUTERMOST INTERSECTIONS OCCJF AT*) 

WRIT£(6*3) 

WRITE! 6*2) CMEGA( I IM I N) , XI C 1 1 Ml N ) * Y I C 1 1 MI N) • 

1X2! J I > , Y2< Jl ) ,CMEGA< J1 ) 

WRITE! 6 *2) OMEGA ( I I MA X) * Xl ( I I M AX ) *Y1 (I I MAX ) * 

1X2! J2>* Y2! J 2 ) • C MEG A! J2 ) 

CHECKING IF F2(W> IS OUTSIDE FI ( W) FOR THE RANGE 
OF OMEGA OF F2 VARYING FROM OKCG A SUB IIMIN TO 
OMEGA SUB UMAX 

WRITE (6*1 ) NF ILE*M * X K , TA Ul • (8CMM ) * MM = 1*5) , (T ( MM > • MM= l *5) 

WRITE (6,40) 

40 FORMAT!/ * 3X * • **C HE CK I NG IF F2 CF OMEGA SJ3 IIMIN THROUGH* * / * 3X, 
1 • F2 OF OMEGA SUB IIMAX LIES OUTSIDE OR INTERSECTS FI***) 

WR ITE ( 6 » 20 ) 

W.R I TE t 6 * 62 ) 

62 FORMAT ( 3X , 

1*L--1 J X2 IS TO THE RIGHT OF FI .UNSTABLE CASE**/*3X* 

2*L=l I F2 INTERSECTS FI, UNSTABLE CASE*./,3X, 

3 * L= 0 : F 2 IS OUTSIDE OF FI, STABLE CASE*) 

WR ITE ( 6 , 55 ) 

55 FORMAT! //, 4X, * OME GA OF F 1 • *6 X* * RE F1*,8X, 

1 * IM FI*, 8X,*R£ F2*,08X**IM F2*, 7X,*0MEGA OF F2 * ,4X * *N • , 3X , * L * ) 

N 3= 1 

N 3M A X— 18-C M*2) 

DO 51 J=IIMIN, IIMAX 

44 CALL SCAN! J * X2 (J)*Y2(J)*X1*Y1 *N , I , l MAX) 

I F ( N ) 46,46,47 

47 L — I. 

52 WR I TE ( 6*53) OMEGA! I > , XI ( I ) * Y1 ( I > ,X2 < J > , Y2 ( J ) , 

40 MEGA ( J ) ,N ,L 

53 FORMAT! 3X * F € • 1 , */2 56 PI * ,3X ,4 !1 PEI 0 . 3 ,3X) ,0PF6 . 1 , '/ 256 PI*,2X,I2, 
1 2 X , 12* * UNSTABLE*) 

GO TO 56 

46 IF( X2( J )-Xl( I ) )71 * 71 ,72 
72 L=-l 

WRITE! 6*53) OMEGA! I) ,Xl ( I ) • Yl (I ) ,X2( J ) • Y2 ( J ) * 

40MEGA! J ) *N *L 
GO TO 56 
71 L = 0 

49 WRITE! 6,48) OMEGA! I) »X1 ( I ) ,Yl II ) ,X2( J)*Y2!J> , 

40MEGA! J ) *N ,L 

48 FORMAT! 3X* F€ • I * * /2 56 PI * *3X ,4 (1 PEI 0*3 *3X) * )°F6.1 , */256 PI • , 2X , I 2, 
l 2 X • 12. * STABLE • ) 

56 IF!N3-N3MAX)54,64, 64 
64 N 3=0 

WR ITE! 6, 1 )NF ILE.M ,XK. TAU1 * ( 8 (MM ) ,MM=i ,5) , (T ! MM ) * MM= 1 ,5 ) 

WR I TE ( 6 , 40 ) 

WRITE! 6, 20) 

WRITE! 6,62) 
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MAIN 


OATE 


7221 * 


16 / 00/ 18 


WRITE! 6.55) 

54 N3=N3*1 
51 CONTINUE 

C PREPARING THE CALCOMP DIAGRAM 

IIIMIN=IIMIN 
ID=IIMAX-I I MI N 
DO 50 1=1. ID 
X1(I) = XU I IMIN) 

YM I >=Y1< I IMIN) 

I IMIN=I IMIN+1 
50 CONTINUE 

CALL SATCHEtXl .Yl. X2 .Y2.ID *JM AX * BUFF ER.^.XK.TAUl t NF1LE. 
1 8 s T . It IMIN. UMAX) 

888 CONTINUE 

CALL PLOT! 0..0. .999) 

999 STOP 
ENO 
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SCAN 


DATE = 7221* 


SUBROUTINE SCANCI • XI .Y1 .X2.Y2 .N. JX * JMAX) 

DIMENSION X2( 1 ) *Y2 ( 1 ) 

SUBROUTINE USED TO FIND POINTS OF INTERSECTION OF 
TWO CURVES. 

DMlN—999. 

DO 3 J=l.JMAX 

D=SQRT( i X1-X2I J) >**2+( Yl-Y2< J > ) * 42 > 

IF(D-DM IN) 2.2.3 

2 DMIN^D 
JX=J 

3 CONTINUE 

IFIDM IN-0.025) € .6. 7 
7 N=0 
RETURN 
6 N=t 
RETURN 
END 


16/08/19 
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SA TCHE 


DATE - 7221* 


16/08/ 10 


SUBROUTINE SA TCHE (XI «Y1 * X2. Y2 • J » It BUFFER* M* XK * TAU 1 * NF ILE * B* 
IT* I IIM IN. I IMAX) 

O IMENSION Xl( 1 ) *Y1 (1 ) *X2( 1) *Y2! 1 >. BUFFER! lhB(H,r(l) 

XL=9. 

YL=7. 

CALL PLOT ( .5. .5*23) 

CALL SCALE! X2.XL.I *1 ) 

X K J + 1 ) = X2 ! 1+1 ) 

Xi( J+2)=X2! 1*2 ) 

CALL SCALE! Y2.YL.1 .1 ) 

Y 1 ! J + l ) — Y2 ! 1 + 1) 

YK J+2)=Y2!I+2> 

C DETERMINING ORIGIN OF AXES FOR X2.Y2 

C ! LOCAT ING HEIGHT OF X-AXIS) 

XSUBX^O. 

IF! Y2! I+l)-0.) 1.2* 2 

2 XSUBY=0.0 
GO TO 5 

1 D=Y2! I+l)+YL*Y2!I+2) 

IF ! D** 0 • ) 3.3.4 
4 XSUBY=-Y2! I+D/Y2! 1+2) 

GO TO 5 

3 X SUB Y= YL 


C LOCATING WIDTH OF Y-AXIS 

5 Y SUB Y= 0 • 

1 F! X 2! I+l)-0.)€.7.7 

7 YSUB X=0 • 

GO TO 10 

6 D=X2! 1 + 1 )+XL*X2! 1 + 2) 
IF!D-0.)8.8.11 

11 YSU8X=-X2( I + l)/X2! 1+2) 

GO TO 10 

8 Y SUB X=XL 

10 CALL AXIS! XSUBX.XSUB Y. • 

1XL* 0. • X2( I + 1) *X2! I +2 ) ) 

CALL AXIS! YSUBX.YSIB Y* • 

1Y2! 1+1 ) • Y2 ! 5+2) ) 

CALL LINE! X2.Y2.I . 1. 0.0) 

CALL LINE! X1.Y1 *J. 1.+3.03) 

I ID = I IMAX- I II MIN 
00 12 JJ=1 . 1 1 D 
X2! JJ) = X2! IIIMIN), 

Y2( JJ)=Y2< IIIMIN) 

I IIMIN=II IMIN+1 

12 CONTINUE 

X2! JJ+1 ) = X 2 ! 1 + 1) 

X2! JJ + 2);=X2! 1+2) 

Y2! JJ + 1 )-Y2! 1 + 1) 

Y2! JJ+2)=Y2! 1+2) 

CALL L INE! X2. Y2 • JJ , 1 .-3 *1 1 ) 
CALL PLOT! YSUBX+O. 20 .YL*23) 
CALL PLOT! 4.8. C.. 22) 

CALL PLOT! 0 ..-1 .1 . 22 ) 

CALL PLOT! — 4 . 8 • 0 . * 22 ) 


REAL F 1 * F2 • .-32. 
IMAG FI *F2* .25.YL.90.. 
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SATCHE 


DATE 


72214 


16/08/ 18 


CALL P LOTI 0.* I •1*221 

CALL SYMBOL! 0 • 15*— 0. 25*0*15 #2 2H S.AT CHE DIAGRAM ,0.,22) 

CALL NUMBER ( l • 6 ,- 0. 4 * O. I • XK • 0. • 2) 

CALL NUMBER (2. 7.-0 .4 .0.1 .TAU1 .0. .2) 

CALL SYMBOL ! 2 . 8 0. 25 . 0 . 1 0 . 1 7H- DEL AY ! S ) P NOBL 5M , 0 . . 1 7 ) 

CALL NUMBER!2. 700, -0.25,0. 10. FLOAT!M) ,0.,-i) 

CALL SYMBOL ! 0. 15.-0. 4.0.1 *3 OHCON ST ANT S : K = .TAU1= ,Q«,30) 

CALL SYMBOL! 0. 15.- C. 55. 0.10 .4 3NT KE EXPONENTS FCP THE PURE TIME DEL 
1 AYS ARE;. 0., 43) 

CALL SYMBOLCO. 15.- 0. 85.0.1 .45HTHE FEEDBACK COEFFICIENTS FOR THE DE 
2LAYS ARE: .0..45) 

XN= 0 . 

DO 1 4 NN=1,5 

C ALLNUMBER! 0. 15 + XN .- 0.7 .0.1 .T ! NN ).0. .1 ) 

XN=XN+0.8 

14 CONTINUE 
XN=0. 

DO 15 NN=l .5 

CALL NUMBER !0.15+XN,-1. .0.1 «8(NN )• 0. • 1 ) 

XN=XN*0.,8 

15 CONTINUE 

CALL SYMBOL!! .2 .-6 .7,0.1 .27HCALC CMP C CMPUT iR-DRAW N PLOT, 0.. 27) 

CALL SYMBOL ! 1 • 2,- 6. 85 ,0.1. 37HPR CGRAM— PLQT DESIGN : JACQUES PRES 

15.0. .37, 

CALL SYMBOL ! I • 2 7. 00 • 0 • 1 , 1 2HDA TE : 8/1 /72 ,0 . • 12 ) 

CALL PLOT! — 4.0,— 0. 2,3) 

CALL PLOT! —3. 7C •— 0 .2 .2 ) 

CALL SYMBOL! ^3.65, -0. 25 . 0 . 1 0 C. 27HF ! EXP! S ) ) : CU RVE WITH DELAYS. 

10., 271 

CALL SYMBOL ! — 4 . 0 .6 5.0. 1*3 OH ++4- + F!S):CURV5 WITH NO DELAYS • 0 ., 30 ) 

CALL SYMBOL(-4 ? .-1 .05,0.1 ,37H**** PART OF F!EXF!S)) CORRESPONDING 

1.0. , 37 ) 

CALL SYMBOLI-3. 65.-1 .20.0.1 .28HTC OMEGA OF F!S) BCUNOEO BY .0..28) 
CALL SYMBOL ! — 3 . 65, — 1 ,35*0,1 ,2 SHI TS OUTER INTERSECTIONS WITH, 0.. 28) 
CALL SYMBOL ! — 3 .65,-1 .50,0.1 ,9HF ! EXPIS)) ,0..9) 

CALL PLOT! 10. .-7.5 .-23) 

RETURN 

END 
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PRECEDING PAGE BLANK NOT EU MED 
APPENDIX F 

THE MATHEMATICAL FUNCTIONS USED IN THE CSMP SIMULATION 

(SEE FIG. 17) 


Computer Notation 

Mathematical Notation 

A. Y = STEP (P) 

Y = 0 t < P 

(STEP FUNCTION) 

Y - 1 t > P 

B. Y = REALP L (IC, P. X) 

PY = Y = X 

Y(0) = IC 
(SIMPLE LAG) 

(nM 

C. Y = INTGRL (IC, X) 

Y = f X dt + IC 

Y (0) = IC 

*a) 

(§) 

(INTEGRAL) 

D. Y = DELAY (N, P, X) 

Y(t) = X(t - P) t ^ P 

P = DELAY TIME 

Y = 0 t < P 

N = NUMBER OF POINTS 


SAMPLED IN INTERVAL P 
(INTEGER CONSTANT) 

( e_sp ) 

DEAD TIME (DELAY) 



Preceding page blank 
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BLOCK 



Figure la. The mathematical box as a component of a system. 
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e - st l = PURE TIME DELAY =TIME LAG 

TS + 1 

Figure 3a. The network graph for a problem with one pure time delay. 


F(s) 
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INPUT 



INPUT 



Figure 4b. The reduction of Figure 4a to two single delays. 


► OUTPUT 


► OUTPUT 
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Figure 6a. The conformal mapping from the s-plane onto the F(s, e s )-plane. 






(A,B are moving points) 


IMAGINARY Fj ,F 2 


INTERSECTION OF 
Fi AND F 2 


F 2 LOCUS 


F(s,e s ) = Fi -F 2 



F 2 (e 5 ): CURVE WITH DELAYS 
— Fj (s): CURVE WITH NO DELAYS 

*(an area is enclosed if it lies to the 
right of the arrow on the locus) 


F7Z: 


RESULTANT VECTOR AB 
AREA ENCLOSED BY F x 


AREA ENCLOSED BY F 2 

Figure 6c. The Sotche diagram. 
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Figure 7. A proposed model for a network flow hindered by two pure time delays. 







>< «r 


Hi=K H2= ^TTT H 3=T D,=e-s«i D 2 = e“ st 2 
INPUT T l s+A 



OUTPUT 

) ► 

*7 


Figure 8- The network graph for the model in Figure 7. 


COMPUTER RUN NUMBER 1 
THE NUMBER OF DECAYS IS 
THE CONSTANTS ARE: 

K* 0 *450 T AU1 = 1 


SATCHE OATA GENERATION 


Bin to e(Mi: 
t c i ) to tcmi: 


i «ooo 

1*000 


1 *000 
2. COO 


OMEGA OF FI 
-256*0/256 PI 
-255*0/256 PI 
-254*0/256 PI 
-253*0/256 PI 
-252*0/256 PI 
-251*0/256 P I 
250*0/256 PI 


RE FI 

-9*87 OE 00 
-9.793 E 00 
-9*71 6E 00 
-9.640E 00 
-9.564E 00 
-9.488E 00 
9.412E 00 


251*0/256 PI 
252*0/256 PI 
253*0/256 PI 
254*0/256 PI 
255.0/256 PI 
256*0/256 PI 


-9.487E CO 
-9.563E 00 
-9.639E 00 
-9.715E 00 
-9.7926 00 
-9.Q69E 00 


NOTATION: 


IM FI 

-3*1 42 E 0 0 
-3*1 29E 00 
-3.1 17E 00 
•3*1 C5E 00 
-3.093E 00 
>3*0e0E 00 
-3.066E 00 


RE F2 
9* C00E-01 
8* 997E-01 
8* 96 6E — 0 1 
6. 97 OE — 0 1 
8* 946E-01 
8. 91 5E-01 



513 LINES ACTUALLY PRINTED 
(COMPRESSED HERE TO SAVE SPACE) 


3.080E 00 
3*0 92 E 00 
3.1 C5E 00 
3*1 17E 00 
3.129E 00 
3*1 41E 00 


8. 91 5E-01 


8* 946E-01 
8. 969E-0 1 
8* 966E-01 
8. 997E-01 
9. 00 OE— 0 l 


RE FI = REAL PART OF F } = X, 

IM FI = IMAGINARY PART OF F = Y, 

RE F2 = REAL PART OF F 2 = X 2 

IM F2 = IMAGINARY PART OF F 2 = Y 2 


IM <=2 

-6 » 26 OE— C 6 
-2 *20 9c — C 2 
-4 .41 5E— 0 2 
-6 .61 7E-0 2 
-3.81 2E-02 
— i 0 Y - 


-0 1 
1 *10 2E-0 1 
8.637E-02 
6 .64 36-0 2 
4.441E-C2 
2.2 35F-02 
2.706E-04 


figure 9. Computer output displaying the generation of numerical 




values for F ^ and F ^ 


OMEGA OF E2 
-256.0/256 ^ I 
-255.0/256 ^ I 
-254.0/256 P I 
-253*0/256 ^ I 
-252.0/256 
-251 .0/256 o l 


. 0 /2 t)6 I 
251 .0/256 » I 
252.0/256 P I 
253.0/256 * I 
2 54.0/256 => I 
2 55.0/256 P I 
256.0/256 ^ I 



COMPUTER RUN NUMBER 
THE NUMBER OF DELAYS 
THE CONSTANTS ARE: 

K« 0.450 TAU1= 


I 

IS 


SATCHE DATA GENERATION 


t .000 


BCD TO e<M>: 
Till TO T (M l: 


1 .000 
I .000 


1.C00 

2.00C 


0 .0 
0 .c 


FINDING ALL POINTS OF INTERSECTION 

OBSERVE THE FOLLOWING CODES BELOW 
N*0 : NOT A REGION OF INTERSECTION 

N*1 : A REGION QF INTERSECTION 


0* 0 

0.0 


OMEGA OF FI 

RE FI 


IM Fl 


RE p 2 

256.0/256 PI 

-9.870E 

00 

-3. I42E 

00 

— 8 * 665 E- 0 1 

255.0/256 PI 

-9.793E 

00 

-3.129E 

00 

-8.665E-0 l 

£56.0/256 PI 

-9.716E 

00 _ 




_ -9.665F-0 I 


7 , r, “ 






(SPACE SAVING) 


0 • 0 

0.0 



OMEf GA OF 
-10.0/256 PI 
-10.0/256 PI 
,0/256 PI 



(ALL ZEROS) 


^0"** . •60.0^'i^ 

ait 

— 5.4 21fc— 0 1 

-7.363E-C1 

— *i • J 1 .. g 1 

— 6 • 37 6 - i. 1 


P 1 

0 

-59.0/256 

pi 

—5 .2426—0 1 

-7.24CF- Cl 

-4.976F-01 

-6. 37 6 £-0 1 

-37 .C/256 

P I 

0 

—58.0/256 

pi 

—5.066c— 01 

-7.118E-CI 

-4.789E-01 

-6.457 £—01 

-3B .0/256 

PI 

0 

-57.0/256 

PI 

-4.693E-01 

-6.CSSE- 01 

-4.789E-0 1 

-6.457 £-0 l 

-38 .0/256 

P I 

0 

-56.0/256 

PI 

-4.723E-01 

-6.672E- Cl 

-4.600E-0 1 

-6.531 £-01 

-39 .0/256 

p I 

0 

-55.0/256 

pi 

— 4.555E— 0 1 

-6.749E- Cl 

-4.410E-C 1 

-6.600 £-0 1 

-40 .0/256 

P 1 

1 

-54.0/256 

PI 

— 4 • 39 IE— 0 1 

— 6. 627E— 01 

— 4.410 E— 0 1 

-6.600 £-01 

-40 .0/256 

P I 

1 

-53.0/256 

p I 

— 4.230E— 01 

-C.5C4E- Cl 

-4.218E-0 1 

-6.631 £-0 1 

-41 .0/256 

p I 

1 

-52.0/256 

pi 

-4.072E-01 

— 6. 38 IE- Cl 

-4.218E-0 1 

-6.661 5-0 1 

-41 iC/256 

PI 

0 

-51.0/256 

PI 

— 3.91 7E— 01 

-6.259E-01 

-4.026E-0 1 

-6.71 7 £-01 

-42.0/256 

PT 

0 

-50.0/256 

pi 

— 3.765E— 0 1 

-6.126E- Cl 

-3.333E-0 1 

-6.766 £-0 1 

—43 .0/256 

P I 

0 

-49.0/256 

pi 

-3.616E-01 

-6.013E-01 

-3.933E-0 1 

-6.766 £-0 1 

-43.0/256 

3 I 

•0 

-48. 0/256 

PI 

-3.670E-01 

-5.890E- 01 

-3.639E-0 1 
— 3.445 “ 

-6.808 £-0 1 

-44.0/256 

P I 

0 


Figure 10. Computer output display ing ‘where Fj and F 2 * n * ersect numerically 


PMFGA OF FI 

PE FI 


P* - 2 

IMF* 

OMEGA OF F 

y 

-♦6.0/256 

PI 

-3.1 375-01 

-5.64SE-C1 

-3.445E-0 1 

-6.844 £-0 1 

-45.0/256 

P I 

-45.0/256 

PI 

-3.049E-01 

-5.522E- Cl 

-3.251E-0 l 

-6.974 £-0 1 

-46 .0/256 

J 1 

-44.0/256 

PI 

-2.9 15E-01 

-5.4C0E- Cl 

-1.3326-01 

-4.750 £-0 1 

-177 .0/256 

P I 

-43.0/256 

PI 

-2.704S-OI 

-3.277E- Cl 

-1 .9326-0 1 

-4.750 £-0 1 

-177.0/256 

P I 

-42.0/256 

PI 

-2.656E-0I 

-5.154E-C1 

-1 .9146-0 1 

-4.619 £-0 1 

-176 .0 /256 

P I 

-41.0/256 

P I 

— 2.5 31 E— 0 1 

-5.031E- Cl 

-1.9146-0 1 

-4.619 £-0 1 

-176. C/256 

P I 

-40.0/256 

PI 

-2.409E-01 

-4.9C9E- Cl 

-l .9146-0 1 

-4.61 9 £-0 1. 

-176 .0/256 

° I 

-39.0/256 

PI 

-2.290E-01 

— 4 • 786E— Cl 

-1.9146-0 1 

— 4 . 61 9 £ -0 1 

-176 .0/256 

® : 

—38 .0/256 

PI 

-2.1 756-01 

-4.663E- Cl 

-1.9896-0 1 

-4.497 £-01 

-175.0/256 

p i 

-37. 0/256 

PI 

-2.062E-01 

-4.«40E- Cl 

-1.9996-0 1 

-4.437 1-0 1 

-175.C/256 

p i 

-36 .0/256 

PI 

— 1 .9526— 0 1 

-4.4ie?-Cl 

-1.9896-0 1 

-4.437 £-0 1 

-175. C/256 

pi 

-35.0/256 

PI 

-l .3456-01 

-4. 2956- 01 

-2.0596-0 1 

-4. 353 £-01 

-174.0/256 

p i 

-34.0/256 

-33.0/256 

V V 

— 1.74 IE— 0 1 
-1 .6406-01 

-4. 172E- Cl 
-4.0506- Cl 

-2.0596-0 1 
-2.123 . 

-4.353 £-0 1 
_ -4. 21 7 £-01 

-174 .0/256 
-1 7 3jj2^j^- 

p i 


(SPACE SAVING) 


(ALL ZEROS) 


<j 

00 


♦ 

-*.0/256 

PI 

-2.405E-03 

-6.135E-C2 
-4.9C8E- C2 

7.149E-0 6 

" v . u 
0.0 

- 1 28 

-126 .0/256 

'pY 

-3.0/256 

PI 

-1 .3E4E-03 

-3.6606- C2 

7.149E-06 

0.0 

-12E .0/256 

P I 

-2.0/256 

PI 

-6.0I8E-04 

-2.453E- C2 

7.149E-0 C 

0.0 

-126 .C/256 

p I 

-1.0/256 

PI 

-I.503E-04 

— 1 • 2 26E- 02 

7.l*9£-<56 

0.0 

-1 26 .0/256 

PI 

0.0/256 

PI 

— 1 .451E— 1 0 

1.205E- 05 

7.149E-06 

0.0 

-1 26 .0/256 


1.0/256 

PI 

— 1 • 509E— 04 

1.228E- 02 

7.149E-C6 

c.o 

-126 .0/256 

P I 

2.0/256 

PI 

— 6.030E— 04 

2.456E- C2 

7.149E-06 

0.0 

-12P .0/256 

P I 

3.0/256 

PI 

—I • 356E-03 

3.663E- C2 

7.1495-0 6 

o 

• 

o 

-126 .0/256 

o j 

4.0/256 

PI 

— 2.41 IE— 0 3 

4.91 0E- C2 

7.149 E— 0 6 

0.0 

-126 .0/256 

P I 

5.0/256 

PI 

-3.7666-03 

6. 127E- C2 

7.149E-06 

0.0 

-126 .0/256 

P I 

6.0/256 

PI 

-5.423E-03 

7.364E- 02 

— 1.101 E— 0 2 

2 . 69 5 £-0 4 

129 .0/256 

P I 

7.0/ 2.56 

PI 

— 7 • 38 IE— 0 3 

e.591E- C2 

— 1.101 E— 0 2 

2.695 £-04 

129 .C/256 

P 1 

8.0/256 

PI 

—9.64 IE— 03 

9.819E- 02 

-2.203E-02 

1 .031 £-03 

1 3C .0/256 

P I 

9.0/256 

PI 

-1 .2206-02 

l.lt5E-Cl 

-2.203E-02 

1.031 £-03 

13C .0/256 

p I 

1 0 .0/256 

PI 

-1 .5C6E-02 

1.227E- 01 

— 3.301 E— 0 2 

2.4J2£-03 

131 .0/256 

p I 

1 \ .0/256 

PI 

-1 .823E-02 

I.3E0E- Cl 

-4.392E-0 2 

4.3 22 £-03 

132.0/256 

P I 

1 2.0/256 

PX 

-2. 169E-02 

1.473F-C1 

-5.474E-02 

* “"**>--0 3 

133 .0/250 

P T 


Figure 10. (Continued) 













COMPUTER RUN NUMBER 1 SATCHE DATA GENERATION 

THE NUMBER OF DELAYS IS 2 
THE CONSTANTS ARE2 
K= 0*450 T AU 1 = 1*000 

ei i ) to b< m ) : !*ooo 1.000 o*o o.o o*o 

Til) TO TIM): 1*000 2*000 0*0 0.0 0.0 


oo THE OUTERMOST INTERSECTIONS OCCUR AT 


OMEGA 3 F F2 
-40.0/256 31 
40.0/256 PI 

Figure 11. Computer output displaying the (numerical values at 
the outermost intersections. 


OMEGA OF FI RE F 1 IM FI R E F2 I M F2 

-55*0/256 PI -4*555E-01 -6*749£-01 -4.410E-01 -6.600E-01 

55*0/256 PI -4.556E-01 6*750E-0l -4.410E-01 6.e03E-01 



COMPUTER RUN NUMBER 2 SATCHE DATA GENERATION 

THE NUMBER OF DELAYS IS 2 
THE CONSTANTS APE: 

K= 0*450 TAU1= 1.000 

B( 1 ) TO B( M ) : l.CCC 1*000 0.0 0*0 0.0 

T( 1 ) TO T ( M ) : 1.0CC 2.000 0.0 0.0 O.C 

★♦CHECKING IF F2 OF OMEGA SUB IIMIN THROUGH 

F2 OF OM^GA SUB UMAX LIES OUTSIDE OR INTERSECTS FI** 

OBSERVE THE FOLLOWING CODES BELOW 
N=0 : NOT A RFGION OF INTERSECTION 

N=1 : A REGION OF INTERSECTION 

L=-l : X2 IS TO THE RIGHT OF FI .UNSTABLE CASE 

L=l : F 2 INTERSECTS FI. UNSTABLE CASE 

L=0 : F2 IS OUTSIDE OF Fl, STABLE CASE 


OMEGA OF 1 

FI 

RE FI 

IM FI 

RE F2 

IM P2 

OMEGA OF 1 

F 2 

N 

L 


-44.0/256 

PI 

-2.915E-C1 

-5.400E-C1 

- 1 .54CE-01 

-6.856F-01 

-55.0/256 

PI 

0 

- 1 

UNSTABLE 

-45.0/256 

PI 

-3.0A9E-01 

— 5 . 522E-C 1 

-1 .724E-01 

-6.883E-01 

-54.0/256 

PI 

0 

-1 

UNSTABLE 

-45. 0/256 

PI 

-3.C49E-01 

-5.522E-C 1 

-l .91CF-C1 

-6.90 3E — C 1 

-53. C/256 

PI 

0 

-1 

UNSTABLE 

-46.0/255 

PI 

-3. 1 87E-01 

-5.645E-C l 

-2.099E-C l 

-6.91 8E-C l 

-52.0/256 

PI 

0 

-1 

UNSTABLE 

-47.0/256 

PI 

-3. 327E-01 

-5.768E-C 1 

-2.2S8E-.01 

-6.926E-C 1 

-51 .0/256 

PI 

0 

-1 

UNSTABLE 

-48.0/256 

R I 

— 3 • 470 E- 0 1 

-5.89CE-C l 

-2.479F-01 

-6.928E-CI 

-50.0/256 

PI 

0 

-1 

unstable 

-48. C/256 

PI 

-3.47CE-01 

-5.890E-C1 

-2.671E-01 

-6.924F-C1 

-49.0/256 

PI 

0 

-1 

UNSTABLE 

-40.0/256 

PI 

-3.616E-01 

-6.013E-C l 

-2.864E-0 1 

-6.91 4E-C1 

-48.0/256 

PI 

0 

-1 

UNSTABLE 

-SC. “5/256 

PI 

-3.765E-01 

-6.136E-01 

-3.057E-G1 

-6.397E-C1 

-47.0/256 

PI 

0 

-1 

unstable 

-50.0/256 

PI 

— 3.7 65 E - 0 1 

—6 • 1 36E-C 1 

-3.251E-01 

— 6.874 F — 0 1 

-46.0/256 

PI 

c 

-l 

unstable 

-5 1.0/256 

PI 

-3.Q17E-01 

-6.259E-C 1 

-3.445E-C 1 

— 6.844E — C 1 

-45.0/256 

PI 

0 

-1 

UNSTABLE 

-52.0/256 

P I 

-4.C72E-01 

-6.381E-01 

-3.639E-C I 

-6.808E-C 1 

-A A. 0/256 

PI 

0 

-1 

UNSTABLE 

-52.0/256 

P I 

— 4.0 72E- 0 l 

—6.381 E- 0 1 

-3.833F-01 

-6.756E-C \ 

-43.0/256 

PI 

0 

-1 

UNSTABLE 

-53.0/256 

PI 

-4.23CE-01 

-6.504E-C1 

-4.C26E-C 1 

-6.71 7E-C 1 

-42.0/256 

PI 

0 

-1 

UNSTAOLT 


Figure 12. Computer output displaying regions of stability and instability 


COMPUTER RUN NUMBER 2 SATCHF DATA GENERATION 

THE NUMBER OF DELAYS IS 2 

THP CONSTANTS ARE: ' 

K= 0.450 TAU1= 1.000 

B( 1 ) TO 9(M): 1.0C0 1.000 0.0 0.0 0.0 

T( 1 > TO T ( M ) I l.CCC 2.000 0.0 0.0 o.c 

♦♦CHECKING IF F2 OF OMEGA SUB IIMIN THROUGH 

F2 OF OMEGA SUB UMAX LIES OUTSIDE OR INTERSECTS FI** 

OBSERVE THE FOLLOWING CODES BELOW 
N=0 : NOT A REGION OF INTERSECTION 

N=1 : A REGION OF INTERSECTION 

L-- 1 : X2 IS TO THE RIGHT OF FI. UNSTABLE CASE 

L=1 : F 2 INTERSECTS Fl.UNSTABLE CASE 

L=0 : F2 IS OUTSIDE OF FI. STABLE CASE 


OMEGA O p FI 

PE FI 

IM FI 

RE F2 

IM F 2 

OMEGA OF 1 

F 2 

N 

L 


-53.0/256 

PI 

— 4 . 230 E — 0 1 

-6.504E-C1 

-4.218E-01 

— 6 * 66 1 F - 0 1 

—4 1.0/256 

PI 

1 

1 

unstable 

-54.0/256 

PI 

-A.391E-01 

-6 .627E-01 

-4.410E-01 

-6.60CE-C I 

-4C. 0/256 

PI 

1 

1 

UNSTABLE 

-55.0/256 

PI 

-4.555E-01 

-6.749E-01 

-4.600E-01 

— 6.531E-C1 

-39.0/256 

PI 

1 

1 

UNSTABLE 

-55.3/256 

PI 

-4.555E-01 

-6.749E-01 

-4.789E-01 

-6.457E-C 1 

-38. C/256 

PI 

0 

c 

STABLE 

-56.0/256 

PI 

-4. 723E-01 

-6.872E-C1 

-4.976E-01 

-6. 376E-01 

-37.0/256 

PI 

0 

r 

STABLE 

-36.0/256 

PI 

-4.723E-01 

-6. 872E-C1 

-5.162E-C1 

— 6 • 26 9E-C 1 

-36.0/256 

PI 

0 

C 

STABLE 

-57.0/256 

PI 

-A.693E-01 

-6.995E-C1 

-5.345E-C 1 

-6. 196E-0 1 

-35.0/256 

PI 

c 

C 

STABLE 

-57.0/256 

PI 

-4.893E-01 

-6.995E-01 

— 5 . 526E-0 1 

-6.096F-C1 

-34.0/256 

PI 

0 

0 

STABLE 

-57.0/256 

PI 

-4.893E-01 

-6.995E-CI 

— 5 . 704E-C; l 

— 5 • 99 1 E— 0 1 

-33. C/256 

PI 

0 

0 

STABLE 

-58.0/256 

PI 

-5.C66E-91 

-7. 1 18E-01 

— 5 • 880E-0 1 

-5.879E-0I 

-32.0/256 

PI 

0 

c 

stable 

-5S. 0/256 

PI 

-5.C66E-01 

-7.1 1BF-01 

— 6 • C52E-0 1 

-E.762E-01 

-31.0/256 

PI 

c 

o 

stable 

-58.0/256 

PI 

-5.C66E-01 

- 7 • 1 1 SE W G 1 

-6.222E-01 

-5.639E-C 1 

-30.0/256 

PI 

0 

c 

STABLE 

-59. 0/256 

P I 

-5. 2A2E-9 1 

-7.2A0E-01 

-6.388E-01 

-5.51 CF-C 1 

-29. C/256 

PI 

0 

0 

STABLE 

-59.0/256 

PI 

-5. 242E-91 

-7.240E-CI 

-6.551E-01 

-5.376E-01 

-28. C/256 

PI 

0 

r 

STABLE 


Figure 12. (Continued). 



COMPUTER run NUMBER 2 SATCHF DATA GENERATION 

THE NUMBER OF DELAYS IS 2 
THF CONSTANTS APE: 

K — 0.450 TAU 1 = 1.0C0 


B( 1 ) TO B( M ) : 1 .OCC 

T ( 1 ) TO T<M)2 l.CCC 


l.OOC O.C 

2. COO 0.0 


0*0 G . 0 

0.0 O.C 


♦ ♦CHECKING IF F 2 OF OMEGA SUB IIMIN THROUGH 

F2 OF OMEGA SUB IIMAX LIES OUTSIDE OR INTERSECTS Fl** 

OBSERVF THE FOLLOWING CODES BELOW 
N=0 : NOT A REGION OF INTERSECTION 

N=1 : A REGION OF INTERSECTION 

L= — 1 : X2 IS to THE RIGHT OF Fl, UNSTABLE CASE 

L=1 : F 2 INTERSECTS F 1 f UNST ABLF CASE 

L = 0 : F 2 IS OUTSIDE OF » ST ABLE CASE 


OMEGA OF ! 

r l 

PF F 1 

IM Fl 

RE F2 

IM F 2 

OMEGA OF 1 

F 2 

N 

l 


-59. 0/256 

PI 

-5. 2425-01 

-7.240E-C1 

-6.709F-01 

-5.2365-0 1 

-27.0/256 

PT 

C 

o 

STABLE 

-60.0/256 

p I 

-5. &21E-01 

-7.363E-C1 

— 6*864E-01 

-5.09 1E-C 1 

— 26.C/256 

PI 

0 

r 

STABLE 

-6 O.C/ 256 

p I 

— 5.421 E -01 

- 7 * 363E- C 1 

-7.01 5E-0 l 

— 4 * 94 C E-C 1 

-25.0/256 

PI 

0 

r 

stable 

— 60. 0/256 

PI 

—5. A21E-01 

- 7 * 363E-C 1 

-7, 161E-C 1 

-4.785E — Cl 

-24.0/256 

PI 

0 

c 

STABLE 

-60.0/256 

p I 

-5.421E-01 

— 7 . 36 3E -0 1 

— 7 • 3C 3E-0 1 

-4.624E-01 

-23.0/256 

PI 

0 

r 

STABLE 

-60.0/256 

p I 

-5*42 1E-0 1 

- 7 * 363E- C 1 

-7.44CE-01 

-4.459F-01 

-22.0/256 

PI 

c 

C 

STA3LE 

-6C.0/256 

p I 

-5.421E-01 

-7. 363E-C 1 

-7.572E-01 

-4 .289E-0 1 

-2 1 .0/256 

PI 

0 

r 

STABLE 

-ec .0/256 

PI 

-5.421E-0I 

-7. 363E-C1 

— 7* 700E-0 1 

-4. 1 1 5E-01 

-20.0/256 

PI 

0 

r. 

STABLE 

-60.0/256 

PI 

-5. 4 2 1 E - 0 1 

— 7 • 363E-C 1 

— 7 , 82 2 E- 0 1 

-3.937E-01 

-19.0/256 

PI 

0 

c 

stable 

-60. 0/256 

PI 

-5.4213-01 

-7.363E-C1 

— 7* 938F- 0 1 

-3.754E-01 

-16.0/256 

PI 

0 

c 

stable 

-60.0/256 

PI 

-5*4215-01 

- 7 ♦ 363F-C 1 

-8.050E-01 

-3.566E-01 

-17. C/256 

PI 

c 

c 

STABLE 

-60.0/256 

p I 

-5*421*3-01 

-7.363E-61 

-8. 155E-01 

-3.378E-C1 

-16.0/256 

PI 

0 

c 

STABLE 

-60*0/256 

p I 

-5.421E-01 

— 7 • 363E -C 1 

— B • 255E-0 1 

-3.184F-01 

-15.0/256 

PI 

0 

c 

stable 

-6C.0/256 

PI 

-5 *42 IE- Cl 

— 7 * 36 3E -C 1 

-8.349E-01 

-2.987E-01 

-14.0/256 

PI 

c 

c 

STABLE 


Figure 12. (Continued). 



COMPUTER PUN NUMBER 2 
THE NUMBER OF DELAYS IS 2 
THE CONSTANTS APE: 

K= C.45C TAU1= 1.0C0 


SATCHE DATA GENERATION 


B(l) TO B < M ) : 1 « C C C 

T ( l ) TO T ( M ) ! l . C C 0 


l.OCC 0* C 

2.000 0.0 


0 . C 0*0, 

0 • 0 C . 0 


♦♦CHECKING IF F2 OF OMEGA SUB IIMIN THROUGH 
F 2 OF OMEGA SUP IIMAX LIES OUTSIDE OR INTERSECTS F 1 ♦* 


OBSERVE THF FOLLOWING CODES BELOW 
N=D ! NOT A PEG I ON OF INTERSECTION 
N=1 : A REGION OF INTERSECTION 

L=-l : X2 IS TO THF r I GHT OF FI. UNSTABLE CASF 

L=1 : F 2 INTERSECTS FI, UNSTABLE CASE 

L=C t F 2 IS OUTSIDE OF FI, STABLE CASE 


OMEGA OF FI 

RE FI 

IM FI 

RE F 2 

IM F 2 

OMEGA OF 1 

F 2 

N 

L 


-6C.C/2S6 

PI 

-5.421E-91 

-7.363E-01 

-8.437E-0 1 

-2.787E-0 1 

-1 3.C/2S6 

PI 

0 

o 

STABLE 

-59.0/256 

PI 

-5. 242E-C1 

-7.240E-C 1 

-8.519E-C 1 

-2.584E-C l 

-1 2. C/256 

PI 

0 

C 

STABLE 

-5.9*0/256 

P I 

-5.242E-01 

-7.240E-CI 

-8.595F-0 1 

— 2.37 8E — C 1 

-1 1 .0/256 

PI 

r 

c 

STABLE 

-59.0/256 

PI 

-5. 2426-01 

-7.24CE-01 

— 3 • 66 5E-0 1 

-2. 1 7CE-0 1 

-10.0/256 

PI 

0 

c 

STABLE 

— 58. C/256 

PI 

-5.C66E-01 

-7. 1 18E-01 

-8.726E-0 l 

-1 .96CE-C1 

-9 • C / 256 

PI 

0 

n 

STABLE 

-53. t /256 

PI 

-5.C66E-01 

-7.1 IQE-C1 

— 8 . 785E- 0 1 

-1.74 7E -G 1 

-8.0/256 

PI 

0 

r 

STABLE 

-57.C/256 

PI 

-4.893E-01 

-6.995E-C 1 

-8.835F-C1 

— 1.53 3E — Cl 

-7. C/256 

PI 

0 

c 

STABLE 

-57. C/256 

P I 

— 4 * £ 9 3E — 0 1 

-6 .995E-01 

— 8 • 679E-0 I 

-1.317E-01 

-6.0/256 

PI 

0 

0 

STABLE 

-56.0/256 

PI 

-4. 723E-C1 

-6.872E-C 1 

-8.916E-0 1 

-1 .C99E-0 1 

-5. C/256 

PI 

0 

c 

STABLE 

-55.0/256 

P I 

-4.555E-0S 

-6.749E-01 

— 3. 946E-0 1 

- 8. 8 0 9 E - C 2 

-4.0/256 

PI 

0 

c 

stable 

-54. C/256 

PI 

-4. 391E-91 

— 6 4 627E-C 1 

— 8 • 970E-0 1 

-6.614E-02 

-3. C/256 

PI 

0 

c 

stable 

-54.0/256 

PI 

- 4 * 39 1 E- 0 1 

-6.627E-C1 

-8.986E-01 

-4.41 3F-02 

-2. C/256 

PI 

0 

0 

stable 

-53.0/256 

P I 

-4.2305-0 I 

-6. 504E-C1 

— 8 • 997E-0 1 

-2.206E-C2 

-1.0/256 

PI 

0 

c 

stable 

52.0/256 

PI 

-4.C72E-01 

6.381E-C1 

— 9# OOOE-O 1 

2. 159E-C5 

0.0/256 

PI 

0 

0 

STABLE 


Figure 12. (Continued). 





COMPUTER PUN NUMBER 2 SATCHE DATA GENERATION 

THE NUMBER OF DECAYS IS 2 
THE CONSTANTS APE: 

K= C.450 TAU1= l.OCO 

B(i) to g(M>: i.oco 1 . 00 c o.o o.o o.c 

T ( 1 ) TO T(M): 1.0CC 2.000 0.0 0.0 c.o 

♦♦CHECKING IF F2 OF OMEGA SUB IIMIN THROUGH 

F2 OF OMEGA SUP UMAX LIES OUTSIDE OR INTERSECTS FI ♦a 

OBSERVE THE FOLLOWING CODES BELOW 
N=0 : NOT A REGION OF INTERSECTION 

N=1 : A REGION OF INTERSECTION 

L = - 1 : X2 IS TO THE RIGHT OF PI .UNSTABLE CASE 

L=1 : F 2 INTERSFCTS FI, UNSTABLE CASE 

L=C : F? IS OUTSIDE OF FI, STABLE CASE 


OMEGA OF FI 

PE FI 

IM FI 

RE F 2 

IM F 2 

OMEGA OF F 2 

N 

L 


53. C/256 

PI 

-4. 23CE-01 

6.504E-01 

-9.997E-0 1 

2.21 IE-02 

1 .0/256 

PI 

0 

c 

STABLE 

54.0/256 

PI 

-4. 392E-01 

6.627E-C l 

-8.986E-01 

4.41 7F-02 

2.0/256 

PI 

C 

c 

STABLE 

54.0/256 

R I 

-4. 3925-0 1 

6.627F-C1 

-8.97CE-0 l 

6.61 BE — C 2 

3.0/256 

PI 

0 

c 

stable 

55.0/256 

PI 

-A.556E-01 

6.75CE-C1 

-8.946E-01 

S.813E-C2 

4.0/256 

PI 

0 

c 

STABLE 

56.0/256 

P I 

-4.723E-01 

6.872E-01 

-8.91 5F-01 

1 . 10CE-C 1 

5.0/256 

PI 

0 

c 

STABLE 

57.0/256 

PI 

— 4.893E- 0 1 

6.995E-01 

-8.878F-CI 

1.317E-C1 

6.0/256 

PI 

c 

0 

stable 

57.0/256 

PI 

-4.893E-01 

6 .995E-01 

-8.635E-01 

1 .533E-0 1 

7.0/256 

PI 

0 

n 

STABLE 

58.0/256 

PI 

— 5. Of 6F-0 1 

7.11 8E — C 1 

-8.7855-01 

1 .746E-C 1 

8.0/256 

PI 

0 

C 

stable 

58.0/256 

PI 

— 5.0 66E“ C 1 

7. 1 1 8E-C 1 

-3.728E-0 1 

1 • 96 0 E — C 1 

9. C/256 

PI 

c 

c 

STABLE 

59.0/256 

PI 

*5. 242E-CI 

7.240E-C1 

-8.665E-01 

2. 171E-C1 

10.0/256 

PI 

c 

c 

STABLE 

59.0/256 

PI 

-5.242E-0 1 

7.24CE“01 

-8.595E-0 1 

2.379E-01 

1 1 .0/256 

PI 

c 

c 

STABLE 

59.0/256 

PI 

-5.2A2C-01 

7.24GE-C1 

-8.519E-C 1 

2.584E-01 

12.0/256 

PI 

0 

c 

STABLE 

60.0/256 

P I 

-5.422F-01 

7.363E-01 

-Q.437E-C1 

2.7B7E-0 t 

l 3.0/256 

PI 

0 

c 

STABLE 

60.0/256 

PI 

-5. A22E“C1 

7.363E-01 

-8.349E-01 

2.988E-C 1 

14. C/256 

PI 

c 

0 

STABLE 


Figure 12. (Continued). 



SATCHF DATA GENERATION 


COMPUTER RUN NUMBER 2 


THF NUMBER OF 

DELAYS IS 

2 




THE CONSTANTS 

are: 





K= 0.450 . 

TAU1= I 

• coo 




B( l > TO B( M) : 

1 .ccc 

1 .000 

0.0 

c.o 

c.o 

T ( 1 ) TO T(M): 

l.OCO 

2.000 

o.c 

0.0 

c.o 


^♦CHECKING IF ^2 OF OMEGA SUB IIMIN THROUGH 

F2 O p OMEGA SUB IIMAX LIES OUTSIDE OR INTERSECTS Fl** 

OBSERVE THE FOLLOWING COOES BELOW 
N=C t NOT A PEG I ON OF INTERSECTION 
N=1 : A REGION OF INTERSECTION 

L=-l : X2 IS TO THE RIGHT OF F 1 .UNSTABLE CASE 

L=1 : F2 INTERSECTS FI* UNST ABLE CASE 

L=C : F 2 IS OUTSIDE OF F 1 . STABLE CASE 

00 


OMEGA 0F Fl 

RF Ft 

IM F 1 

RE F 2 

IM F 2 

OMEGA OF F2 

N 

L 


60. C/256 

PI 

-5.422E-01 

7. 363E-C 1 

-8.255E-01 

3. IS5E-C1 

15.0/256 

PI 

0 

c 

STABLE 

60.0/256 

P I 

-E.422E-01 

7.363E-C1 

-8.155E-01 

3.378E-01 

16.0/256 

PI 

C 

c 

STABLE 

6C .0/256 

PI 

— 5 • A 22E-0 1 

7.363E-C l 

— 8 * 049E-C l 

3.566E-C 1 

17.0/256 

PI 

0 

c 

stable 

6C. C/256 

PI 

-5.422E-01 

7.363E-C1 

-7.938E-01 

3.755E-0 1 

18. C/256 

PI 

0 

c 

stable 

60.0/256 

PI 

•5. A22E-0I 

7.363E-C1 

-7.821E-01 

3.937E-01 

19.0/ 256 

PI 

0 

i"* 

STABLE 

6C. 0/256 

PI 

-5.422E-01 

7.363E-C 1 

-7.699E-01 

4.U6E-01 

20. C/256 

PI 

0 

r 

stable 

60.0/256 

PI 

-5.422F-01 

7.363F-01 

- 7 . 572F-Q 1 

4 . 290 E — C 1 

2 1.0/256 

PI 

0 1 

c 

STABL E 

60.0/256 

PI 

-5. 4225-01 

7.363E-01 

— 7 . 440E-0 1 

4.4E9E-01 

22.0/256 

PI 

0 

0 

STABLE 

60. 0/256 

PI 

-5. A?2E-01 

7.363E-C1 

— 7 # 30 3E-C 1 

A.625E-C1 

23.0/256 

PI 

0 

c 

stable 

60*0/256 

PI 

-5.422E-01 

7.363F-01 

— 7.161 E — 0 1 

4.785E-01 

24 .0/256 

PI 

c 

r\ 

stable 

60.0/256 

PI 

-5.422E-01 

7 • 363F-G 1 

-7.015E-01 

4.94CE-01 

25.0/256 

PI 

0 

r 

stable 

60 .0/256 

PI 

-E. A22E-01 

7.363E-C1 

— 6 • 864E— C 1 

5.09 IE-0 1 

26.0/256 

PI 

c 

c 

STABLE 

59.0/256 

PI 

«5.242E-0 l 

7 . 240F-0 1 

— 6 • 7C9E-0 1 

5.236F-C1 

27. C/256 

PI 

0 

c 

stablf 

5 9. C/ 256 

PI 

-5.2A2E-01 

7. 240E-C 1 

-6.550E-C 1 

5.376E— 01 

2 8.0/256 

PI 

0 

r 

STABLE 


Figure 12. (Continued). 


COMPUTER RUN NUMBER 2 

the number df delays is 2 

THE CONSTANTS arf: 

K= 0.450 TAU1= 1.000 


SATCHF DATA GENERATION 


8(1) TO 8( M ) : 
T ( 1 ) TO T ( M ) : 


i . f c c 
l .ccc 


l • c o o 

2.000 


♦ ^CHECKING IF F P OF OMEGA SU0 IIMIN THP OUGH 

F2 OF OMEGA SUP UMAX LIES OUTSIDE OR INTERSECTS Fl** 

OBSERVE THE FOLLOWING CODES BELOW 
N= ? : NOT A REGION of INTERSECTION 

N=1 : A REGION OF INTERSECTION 

L=— 1 : X2 IS TO THE RIGHT OF FI. UNSTABLE CASE 

L=1 : F 2 INTERSECTS FI, UNSTABLE CASE 

L=C : F 2 IS OUTSIDE OF Fl , ST ABLE CASE 


3MEGA OF 

Fl 

RE Fl 

IM Fl 

PE F2 

IM F? 

OMEGA OF 

F2 

N 

L 


59.0/256 

D I 

-5.2A2E-C1 

7. 24CE-0 1 

-6.388E-01 

5.510E-01 

29.0/256 

PI 

0 

c 

STABLE 

88. 0/256 

PI 

-5.C 66E-01 

7.1 18E-C1 

-6.222F-01 

5.639E-C 1 

30 .C/256 

PI 

C 

r 

STABLE 

58.0/256 

PI 

— 5 . C 66E — : I 

7. 1 18E-C1 

-6.C52E-C1 

5.762E-0 1 

31 .0/256 

PI 

c 

n 

STABLE 

58.0/256 

PI 

-5.C66E-01 

7. 1 18F-01 

-5.879E-01 

5.8SCE-C 1 

32.0/256 

PI 

c 

r. 

STABLE 

57.0/256 

PI 

-4.893E-01 

6.995E-C1 

-5.704E-01 

5.99 IE — Cl 

33.0/256 

P I 

0 

C 

STABLE 

57.0/256 

PI 

-A. P93E-01 

e .995E-C 1 

— 5 . 525E-T 1 

6.097E-C1 

34.0/256 

PI 

0 

c 

stable 

57. C/256 

PI 

— A . 99 3E- 0 I 

6.995E-C1 

-5. 344E-0 1 

6. 196E-01 

35.0/256 

PI 

0 

r* 

STABLE 

56.0/256 

PI 

-4.723E-01 

6.872E-01 

-5. I61E-01 

6.289E-C1 

36.0/256 

PI 

0 

c 

STABLE 

56.0/256 

PI 

-4.723E-C1 

6 .872E-C 1 

-4.976E-01 

6.376F-C 1 

37.0/256 

PI 

c 

C 

STABLE 

55.0/256 

P I 

-4. 556E-01 

6.750E-01 

— 4 * 78 9F — 0 1 

6.457E-C1 

38.0/256 

PI 

0 

r 

STABLE 

55. C/256 

p I 

-4. 556E-01 

6.750E-C1 

-4.60CE-C 1 

6.5 32E — 0 1 

39.0/256 

PI 

1 

l 

UNSTABLE 

54.0/256 

P I 

-4. 392E-01 

6.627F-01 

-4.410E-C1 

6.600E-C1 

40.0/256 

PI 

1 

l 

UNSTABLE 

53.0/256 

P I 

-4. 2306-01 

6.504E-0I 

-4.218E-01 

6.66 1E-0 1 

41.0/256 

PI 

1 

l 

UNSTABLE 

53. C/256 

P I 

-4. 230 E-01 

6.504E-01 

—4 • C26F-0 1 

6.717E-C 1 

42.0/256 

PI 

0 

- 1 

UNSTABLE 


Figure 12. (Continued). 



SATCHE DATA GENERATION 


00 

00 


COMPUTER RUN NUMBER 2 
THE NUMBER OF DELAYS IS 2 
THE CONSTANTS ARF : 

K= 0.450 T AU 1 = 1.000 

B( 1 ) TO B < M ) : l.CCC 1 .000 0.0 0.0 0.0 

t( n to tim): i.coc 2.000 0.0 0.0 c.o 

♦♦CHECKING IF F2 OF OMEGA SUB IIMIN THROUGH 

F2 OF OMEGA SUB UMAX LIES OUTSIDE OR INTERSECTS FI** 

OBSERVE THE FOLLOWING COOES BELOW 
N=0 : NOT A REGION OF INTERSECTION 

N=1 : A REGION OF INTERSECTION 

L=- 1 : X 2 IS TO THE RIGHT OF FI. UNSTABLE CASE 

L=1 : F 2 INTFRSFCTS FI, UNSTABLE CASE 

L=0 : F 2 IS OUTSIDE OF FI, STABLE CASE 


OMEGA OF FI PE FI 


52.0/256 

P I 

-4.C72E-01 

52.0/256 

PI 

— 4.C72E-01 

5 1 . C/256 

P I 

-3.917S-01 

50. C/256 

PI 

-3. 765E-01 

50.0/256 

PI 

-3.765E-01 

49.0/256 

P I 

-3.616E-C1 

48.0/256 

PI 

-3.47CE— 01 

46.0/256 

PI 

-3. 470E-0 1 

47.0/256 

PI 

— 3. 3 2 7 F — 0 1 

*6.0/256 

PI 

-3. 187E-01 

*5.0/256 

PI 

— 3 • 0 50E - 0 1 

45.0/256 

PI 

-3. C5CE-01 

44.0/255 

PI 

-2.916E-01 


IM FI 

RE F2 

6.381 E— 0 1 

-3.832E-C1 

6.381E-01 

-3.639E-0 1 

6.259E-01 

-3.445E-01 

6 . 1 36E-01 

-3.251E-C1 

6.1 36E-01 

-3.C57E-01 

6.C13E-C1 

-2.864E-01 

5.391E-C1 

-2.671E-01 

5.891E-C 1 

-2.479E-01 

5.768E-CI 

-2.288E-C 1 

5.645E-C 1 

-2.098E-C1 

5.522E-C 1 

-1.91 0E-C 1 

5.522E-C1 

- 1 .724E-0 1 

5.40CE-01 

— 1 • 539E-0 1 


I M F 2 
6 .7S6E-C 1 
6.3C8E-G1 
6.8 4 4E — C 1 
6.S74E-01 
6.S9TF-C 1 
6.91 4E-C 1 
6.924E-C 1 
6.928E-C 1 
6.926E-C 1 
6.918E-C1 
6 • 9 0 3E -0 1 
6.33 2E-C 1 
6.656E-C1 


OMEGA OF 1 
43.0/256 
44.0/256 
45.0/256 
*6.0/256 
47.0/256 
46.0/256 
49.0/256 
50.0/256 
51.0/256 
52.0/256 
53.0/256 
54.0/256 
55.0/256 


2 N L 

PI 0-1 

PI 0 - 1 

PI 0-1 

PI 0-1 

PI 0 - 1 

PI 0-1 

PI 0-1 

PI 0 -1 

PI 0-1 

PI 0-1 

PI 0-1 

PI 0-1 

PI 0-1 


UNSTABLE 

UNSTABLE 

UNSTABLE 

UNSTABL 

unstabl 

UNSTABLE 

UNSTABLE 

UNSTABLE 

UNSTABLE 

UNSTABLE 

UNSTABLE” 

UNSTABLE 

UNSTABLE 


Figure 12. (Concluded) 
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O.UO 0.8d 1.20 

REAL FJ,F2 


CALCQMP COMPUTER-ORAMN PLOT 
PROGRAM-PLOT DESIGN i JACQUES PRESS 
OATEi 8/1/72 


Figure 13. First Satche diagram for a case with two delays. 
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o 


OJ 



Figure 14. Second Satche diagram for a case with two delays. 
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-FlEXPtSn iCURVE WITH DELAYS 


♦♦♦♦ If (S) (CURVE WITH NO DELAYS 

o 

HHHH PART OF FJEXPtSH CORRESPONDING^- 
TO OMEGR OF F (S) BOUNOED BY 
ITS OUTER INTERSECTIONS WITH 
FCEXP(S)) 



.20 


Figure 15. Third Satche diagram for a case with two delays. 



o 

® 



Figure 16. Fourth Satche diagram for a case with two delays. 
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****tONTINUOl,S SYSTfW vnr)ELl N'*, ppCGPAM**** 


*** VERSION 1 .3 *** 

LABEL TIME O^LAY MODEL UNIT ST C F INPUT 
CONSTANT 81 = 1 . .02=1. 

PARAMETER XK = 0 • 4 5 * TAU1=1 

T 1 = 1 • ♦ T2 = 2. 

D YN AM IS 
X1=ST«(0. ) 

X2=X1-B 1* X6-R?*X7 
X3=XK*X? 

X4=PEA^^L ( o. • taui . XI > 

X«= IN TGRL (0«|X4) 

XE=DELAY( 1 0# T 1 * X5I 
X^=DELA Y( 10.T2i X6 ) 

T IMFP ^EL T=0 • 05t E I N TIM = P0* *PR">PL =1 . 0 *QUTr)FL = 1. 0 
METHOD ADAMS 
PRTPL T X7 
ckjh 

P apamftfp xk-o.25 (PROGRAM IS RERUN AT THIS STEP) 

END 

parameter xk = o • 3 5 (PROGRAM IS RERUN AT THIS STEP) 

END 

STOP 

OUTPUT VAPTARLE SEQUENCE 

X 6 X7 XI X? X 3 7^7 0 0 03 XA X5 

OUTPUTS TNP'JTS p AP AMS INTPGS 4- MEV RLKS FORTRAN DATA C )S 

12(500 ) 3R < 1400) 1 0( 400) ?♦ 2= 4( 300) 3(500 12 


END JG R 


NOTATION: 

XK= K (AS IN TEXT) 

REALPL, INTGRL, DELAY, STEP = FUNCTIONS DESCRIBED IN APPENDIX F. 

METHOD ADAMS = INTEGRATION METHOD USED. (SEE REF. 52) 

TIMER = DESCRIBES THE SELECTED TIMING SEQUENCE. 

Figure 17. Continuous System Modeling Program (simulation) for the case of two delays. 
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TIME 0ELA7 Ml DEL UNIT STEP INPUT 


TIME 


XT 

0 .0 


0.0 

1 • OOOOE 

00 

0.0 

2. OOOOE 

00 

0.0 

3 * OOOOE 

00 

o.o 

A .OOOOE 

00 

1 .6572E-01 

5. OOOOE 

00 

5.0530E-01 

6 .OOOOE 

00 

8.5712E-01 

7. OOOOE 

00 

1.112BE 00 

8. OOOOE 

00 

1.1848E 00 

9 .OOOOE 

00 

1.0273E 00 

1 .OOOOE 

01 

6. 745- SE- 01 

1 . 1000E 

01 

2 . 3200E— 0 1 

1 • 2000E 

01 

- 1.5673E-01 

1 .3000 E 

01 

- 3. 5368E-01 

1 • 4000 E 

01 

- 2. 7464E-01 

I . 50 00E 

01 

7.701 4E- 02 

1 .6000E 

01 

6. 0409E-01 

1 .7000 E 

01 

1.1408E 00 

1 .8000 E 

01 

1.5039E 00 

I .9 000 E 

01 

1.5 5 3 BE 00 

2 .OOOOE 

01 

1.2464E 00 

2*1 000 E 

01 

fi • 56 37E- 0 1 

2 • 2000 E 

01 

- 3.9105E-02 

2.3000E 

Cl 

-6.0776E-01 

2 * 4000 E 

01 

-0. 4637F-01 

2 .5000E 

01 

-6. 4355E-01 

? * 60 00 E 

01 

-3.0873E-02 

2 • 7 000 E 

01 

8. 1673E-01 

2 .8000 E 

01 

1.6286E 00 

2 *9000 E 

01 

2.1232E 00 

3. OOOOE 

01 

2.1034E 00 

3.1 000 E 

01 

1.5317E 00 

3.2000E 

01 

5 • 56 T 6E— 0 1 

3 • 3000F 

01 

- 5. 2163E-01 

3 *4 000 E 

Cl 

- 1 • 3424E 00 

3 • 5000 E 

01 

-1.6O10E 00 

3 .6000E 

01 

-1.1620E 00 

3 .7000E 

01 

- 1 . 1433F- 01 

3 .8000 E 

01 

1.2344E 00 

3 .9000 E 

01 

2.4473E 00 

4. OOOOE 

01 

3.0965E 00 

4 . 1000E 

01 

2.9106E 00 

4 .2000E 

01 

1.P816E 00 

4 • 30 00 E 

01 

2. 8990E-01 

4 *4000 E 

01 

-1.3669F 00 

4 .5000 E 

01 

-2.5256E 00 

4 .60 00 E 

01 

-2.7461E 00 

4 .7000E 

01 

-1.0717E 00 

4 .POOOE 

01 

- 1 .C819E-01 

4 .9000 E 

01 

2.01 5°F 00 

5 .OOOOE 

01 

3.8O20E 00 


MINIMUM 
-1.2101E 01 
I 


VERSUS TIM* 








Figure 18. Output to Figure 17 with K -0,45- 
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TIME DELAY 


TIME 

5.1000E 01 
5.2000E 01 
5 . 3000 E 01 
5.4000E 01 
5.5000E 01 
5.6000E 01 
5.7000E 01 

5 *3000 E 01 
5 • 90 00 E 01 
6.0000E 01 
6.1000E 01 
6.2000E 01 
6.3000E O! 
6 • 4000 E 01 

6 « 5000 E 01 
6.6000E 01 
6.T000E 01 
6.8000E 01 
6.9000E 01 
T.OOOOE 01 
7.1000E 01 
7.2003E 01 
7 . 3000 E 01 
7 • 4000 E 01 
7.5000E 01 
7.6000E 01 

7 *70 00 E 01 
7.8000E 01 
7.9000F 01 
9.0000E 01 


MODEL UNIT 


X7 

4.61 03E 
4.076QE 
2.2703E 
- 2.9934E- 

- 2. 3160E 

- 4 • 4 1 1 7E 
-4.4S94E 

- 2. 8095E 
1.181 1E- 
3. 4289E 
6. 0198E 
6.9399E 
5.7295F 
2.6204E 

- 1 . 4840E 

- 5. 261 IE 
-7.3868E 

- 6.9893E 
-3.9877P 

8. 1 213E- 
5.91 84E 
9 .6052E 

1 . 0485E 
^8 *-0X4 3 E 
' 2 • 753 IE 
- 3.7349E 

- 9. 3286E 

- 1 • 2032E 

- 1 .0667E 
-5.3S38E 


STEP INPUT 

MINIMUM 
-1.2101E 01 
1 

00 

00 

00 

01 

00 

00 

00 

00 

01 

00 

00 

00 

00 -- 

00 

00 

00 

00 

00 

00 

01 

00 

00 

01 

00 

00 

00 


00 ♦ 

01 + 

01 + 

00 



18. (Continued) 
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TIME DELAY M3 DEL UNIT STEP INPUT 


TIME 


X7 

0 .0 


0.0 

1 .0000 E 

00 

0.0 

2 .0000E 

00 

0.0 

3 • OOOOE 

00 

0.0 

4 .OOOOE 

00 

9.2065E-02 

5. OOOOE 

00 

2.821 4E-01 

6 .OOOOE 

00 

4.9230E-01 

7. OOOOE 

00 

6.7932S-01 

9 .OOOOE 

00 

8. 0689E-01 

9 .OOOOE 

00 

8.5620E-01 

1 .OOOOE 

01 

8. 2467P-01 

1 . 1000E 

01 

7.2930E-01 

1 . 200 0 E 

01 

5. 9470E-01 

1 . 3000 E 

01 

4. 5623E-01 

1 . 40 00 E 

01 

3.4264E- 01 

I .5000 E 

01 

2.7509E-01 

l .60 00 E 

01 

2.6257E-01 

1 . 7000 E 

01 

3.0153E-01 

1 .8000 E 

01 

3.7794E-01 

1 .9000E 

01 

4.7144E-01 

2 .OOOOE 

01 

5.6039E-01 

2 . 1000E 

01 

6. 2655E-01 

2 . 2000 E 

01 

6.5355E-01 

2 . 3000 E 

01 

6.S350E-01 

2 .4000 E 

0! 

6.1666E-01 

2 . 5000 E 

Cl 

5. 5939E-01 

2 . OOOOE 

01 

4. 9607E-01 

2.7000E 

01 

4.4074E-01 

2 • 8000E 

01 

4. 041 6E-01 

2 .9000 E 

01 

3.9192E-01 

3. OOOOE 

01 

4.0385E-01 

3 . IOOOE 

01 

4. 34 70E-01 

3 • 2000 E 

01 

4.7584E-01 

3 . 3000 E 

ot 

5. 1741E-01 

3.4000E 

01 

5.5062E-01 

3 . 5000 E 

01 

5.6940E-01 

3 . 6000 E 

01 

5.7142E-01 

3 .7000E 

01 

5. 5816E-01 

3.9000E 

01 

5. 3413E-01 

3 .9000E 

01 

5.0563E-01 

4 .OOOOE 

01 

4.7916E-01 

4. IOOOE 

01 

4.6007E-01 

4 . 200 0 F 

01 

4.515^5-01 

4.3000E 

01 

4.5423E-01 

4 • 4000E 

01 

4.6625E-01 

4 . 5000 E 

ot 

4 * P402E— 01 

4 . 6000 E 

01 

5.031 5E- 01 

4 • 7Q00E 

01 

5. 1945E-01 

4 . 8000 E 

01 

5.2981 E— 0 1 

4 .9000 E 

01 

5.3271 E— 0 1 

5 .OOOOE 

01 

5. 2339E-01 


MINIMUM 

0.0 

I 


VERSUS TIME 


MAX IMUM 
3 .5620E-0 1 
I 








+ 

— — -+ 


Figure 19. Output to Figure 17 with K = 0.25. 
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TIME DELAY M3 DEL UNIT STEP INPUT 


TIME 

5.1000E 01 
5 • 20 00 E 01 
5 • 30 00 E 01 
5.4000E 01 
5 ■ 5000 E 01 

5 *60 00 E 01 
5.7000E 01 
5.8000E 01 
5 • 9000 E 01 
6.0000F 01 
6.1000E 01 
6 « 2000 E 01 

6 * 3000 E 01 
6 • 4000 E 01 
6.5000E 01 
6 • 600 0 E 01 

6 .7000E 01 
6.8000E 01 
6.9000F 01 
7 • 0 000 E 01 

7 « 1000 E 01 
7 • 2000 F 01 
7 • 30 00 E 01 
7.4000E 01 
7.5000E 01 
7 • 60 00 E 01 
7.7000E 01 
7.8000E 01 
7.9000E 01 
3 .0000E 01 


X7 

5. I 857E-01 
5, 0594E-01 
4. 9349E-01 
4 « 8 382E- 01 
4.7965E-01 
4.7857E-01 
4. 830 IE-01 
4.9054E-01 
4. 9920E- 01 
5*07 06E— 0 1 
5 • 1 254E- 01 
5 • 1475E-01 
5* 1 361E- 01 
5 • 0972F— 0 1 
5. 0423E-01 
4.9846E-01 
4 • 9 367E— 0 1 
4 • 9076E— 0 1 
4,901 2E-01 
4, 91 64E-0 1 
4.9475E- 01 
4. 9861F- 01 
5,023 3E— 0 1 
5,051 4E-01 
5.0655E- 01 
5,064 IE— 0 1 
5, 0494E-01 
5, 0259E-01 
4.9997E-01 
4.9764E-01 


MINIMUM 
0 , 0 

I 


VERSUS TIM; 


Figure 19. (Continued) 
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DELAY M3 DEL UNIT STEP INPUT 


TIME 

0 .0 

1 • 0 0 OOE 00 
2.0000E 00 
3.0000E 00 
4.0000E 00 
5.0000E 00 
6 • 0 00 0 E 00 
7 • 00 00 E 00 
9.0000F 00 
9 • 0 0 00 E 00 
1 • 00 00 E 01 
1 . 1000E 01 
1.2000E 01 
1 . 30 00 E 01 
1 .4000E 01 
1 • 50 00 E 01 
1 *60 00 E 01 
1 • 7000 E 01 
1 • BO 00 E 01 

1 • 9 0 00 E 01 
2 . 00 00 E 01 

2 • 1 0 0 0 E 01 
2.2000F 01 
2 • 30 00 E 01 
2* « 4000 E 01 
2 • 50 00 E 01 
2 • 60 OOE 01 
2.7000E 01 
2 • 8000E 01 

2 *9000E 01 

3.0000E 01 
3.1000E 01 
3 • 20 00 E 01 
3.3000E 01 
3.4000E 01 

3.5000E 01 
3 • 60 00 E 01 
3.7000E 01 
3.R000E 01 

3 * 9000 F 01 

4 .0000E 01 
4.1000E 01 
4.2000E 01 
4.3000E 01 
4 • 4000 E 01 
4 • 50 OOE 01 
4 • 6000 E 01 
4 • 7 000 E 01 
4.8000E 01 
4.9000E 01 
5.0000E 01 


MINIMUM 
•2. 15 31 E— 02 
I 


VERSUS TIME 


M A X I MUM 
1 .0 30 6E 0 0 
I 


Figure 20. Output to Figure 17 with K =0,35- 
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TIME DELAY MD DEL UNIT STEP INPUT 


T IME 

5*1 000 E 01 
5 • 200 0 E 01 
5.3000E 01 
5.A000E 01 
5.5000E 01 
5 .6000E 01 
5 • 7 000 E 01 
5.8000E 01 
5.9000E 01 
6.0000E 01 
6.1000E 01 
ft • 20 0 0 E 01 
6.3000E 01 
6 • 4000 E 01 
6 • 5000 E 01 

6 .6000E 01 
6.7000E 01 
6.8000E 01 
6.9000E 01 
7.0000E 01 
7,1 0 00 E 01 
7,20 00 E 01 
7 , 30 00 E 01 
7 , 40 00 E 01 
7,'SOOOE 01 
7.60 00 E 01 

7 * 7000 E 01 
7 , 8 000 E 01 
7 ,90 00 E 01 
3 ,0000 E 01 


UNIMUM 
, 1 531E-02 
T 


X7 VERSUS TI Mi 


MAXIM UM 
1.C306E 00 
I 
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Figure 21. First Satche diagram for a single delay case. 
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Figure 22. Second Satche diagram for a single delay case. 
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Figure 23. Third Satche diagram for a single delay case. 
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f^lEXP (S) ) * CURVE WITH DELAYS 

♦ ♦♦♦ (S) i CURVE WITH NO DELAYS 

mw PART OP FJEXPIS)) CORRESPONDING 
TO OMEGA OF F (S) BOUNDED BY 
ITS OUTER INTERSECTIONS WITH 
F (EXP IS)) 



Figure 24. Fourth Satche diagram for a single delay case. 
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Figure 25. Fifth Satche diagram for a single delay case. 
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Figure 26 . First Satche diagram for a case with three delays. 
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Figure 27. Second Sotche diagram for a case with three delays. 
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Figure 28 . Third Satche diagram for a case with three delays. 
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Figure 29. Satche diagram for a case with four delays. 
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Figure 30. First Satche diagram for a case with five delays. 
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Figure 31. Second Satche diagram for a case with five delays. 


110 




Ill 


=p 



Figure 32. Third Satche diagram for a case with five delays. 




